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ABSTRACT 

A noise-based non-parametric technique for detecting nebulous objects, for example, irregular or clumpy galaxies, 
and their structure in noise is introduced. “Noise-based” and “non-parametric” imply that this technique imposes 
negligible constraints on the properties of the targets and that it employs no regression analysis or fittings. The sub¬ 
sky detection threshold is defined and initial detections are found, independently of the sky value. False detections 
are then estimated and removed using the ambient noise as a reference. This results in a purity level of 0.89 
for the final detections as compared to 0.29 for S Extractor when a completeness of 1 is desired for a sample of 
extremely faint and diffuse mock galaxy profiles. The difference in the mean of the undetected pixels with the known 
background of mock images is decreased by 4.6 times depending on the diffuseness of the test profiles, quantifying 
the success in their detection. A non-parametric approach to defining substructure over a detected region is also 
introduced. NoiseChisel is our software implementation of this new technique. Contrary to the existing signal-based 
approach to detection, in its various implementations, signal-related parameters such as the image point spread 
function or known object shapes and models are irrelevant here. Such features make this technique very useful in 
astrophysical applications such as detection, photometry, or morphological analysis of nebulous objects buried in 
noise, for example, galaxies that do not generically have a known shape when imaged. 

Key words: galaxies: irregular - galaxies: photometry - galaxies: structure - methods: data analysis - techniques: 
image processing - techniques: photometric 

Free reproduction: All the data-generated numbers and figures in this paper are exactly reproducible/configurable 
with a make command. See reproduce /readme in the arXiv source files. All results generated by free software. 


1. INTRODUCTION 

Galaxies are one of the most prominent forms of nebulous 
or amorphous objects in astrophysics and astronomical image 
processing. They appear to have a very diverse kinematic his¬ 
tory and thus display a very large variety of shapes and forms. 
They also have complex three-dimensional structures which we 
can observe in two dimensions. Galaxy images can host an ar¬ 
bitrary number of clumps positioned anywhere over their light 
profiles. These clumps might be minor or major mergers, su¬ 
pernovae, other galaxies, or stars on the same line of sight, or 
highly localized unobscured star-forming regions. With such 
hurdles, astronomical image processing, and in particular galaxy 
morphological analysis, can be considered one of the most de¬ 
manding of all disciplines in image processing. 

Observations are inevitably diluted with noise. Any system¬ 
atic bias in measuring the noise properties will directly propagate 
to all higher level measurements on the detected signal or scien¬ 
tific targets. Therefore, accurately measuring the noise charac¬ 
teristics is the most fundamental and primary issue in data analy¬ 
sis and subsequent astrophysical interpretations. However, noise 
can only be accurately characterized when the signal or object(s) 
buried in the noise are detected and removed. Therefore any dif¬ 
ficulties in detection can directly translate into a systematic bias 
in the measured noise properties. For example, objects that are 
very faint and diffuse overall are very hard to detect. For brighter 
objects that can be detected, identifying the possible boundaries 
is a major concern. The uniformity in shape, or morphology, 
between various objects can also play a significant role in the 
detection ability. 

Detection of nebulous astronomical targets significantly suf¬ 
fers from all the problems mentioned above. Because of the in¬ 


strument and atmospheric point spread function (PSF) and also 
the intrinsic shapes of galaxies, their profiles sink deep into the 
noise very slowly with a hardly detectable clear cutoff. Astro¬ 
nomical targets also have an extremely wide range of appar¬ 
ent luminosity and surface brightness profiles, from very bright 
nearby stars and galaxies to faint objects far below the detection 
limits. 

The most commonly used detection methods in astronomy, 
regardless of their implementation, can be classified as signal- 
based detection (see Appendix B for a review). The signal’s a 
priori known properties are the basis of this method. Thus the 
inherent shape (for example, being an ellipse or a functional ra¬ 
dial profile) of the galaxies and the atmospheric and instrumental 
effects, or PSF, have to be known prior to running the detection 
algorithm. This approach to detection will therefore be most ac¬ 
curate for objects that satisfy the a priori known shapes and will 
become less reliable as the inherent shapes of the targets deviate 
from the expected shape. Failing to correctly find the PSF will 
also hamper the accuracy of the results. 

Systematic biases that are caused by cosmological surface 
brightness dimming along with K corrections and morphologi¬ 
cal evolution can hamper the accuracy of any study attempting 
to compare galaxy populations at different redshifts. Therefore it 
is very important that the detection and measurement tools used 
can provide less biased results. An example of such compar¬ 
isons is the observed morphological evolution of blue galaxies 
where it has been observed that at higher redshifts (z > 1), a 
larger fraction of galaxies display irregular or clumpy morpholo¬ 
gies compared to their lower-z counterparts (Cowie et al. 1995; 
Abraham et al. 1996; Lotz et al. 2006; Elmegreen et al. 2007; 
Murata et al. 2014). This observed morphological evolution bi¬ 
ases the photometric measurements of any detection method that 
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is based on a fixed elliptical shape and surface brightness pro¬ 
file for all redshifts, for example, the Petrosian (1976) and Kron 
(1980) methods. It is currently thought that the the peak of cos¬ 
mic star formation density was at z ^ 2 (Hopkins and Beacom 
2006; Kajisawa et al. 2010; Lundgren et al. 2014). Hence, the 
spatial variation of star formation activity within these earlier 
galaxies or their star formation morphology is instrumental in 
understanding star formation in the galaxies. 

The concordance or ACDM cosmological model is another 
field that heavily depends on detection ability. It predicts much 
fainter, lower surface brightness dwarf galaxies than currently 
observed and may be lurking below the currently used detection 
thresholds (Kauffmann et al. 1993; Klypin et al. 1999). From the 
theoretical point of view, changing the astrophysical constraints 
of the semi-analytic models results in the simulations and ob¬ 
servations becoming more consistent (see Maccio et al. 2010, 
and references therein). However, an observational confirma¬ 
tion is still lacking. Since such low mass galaxies are very faint 
and have irregular/amorphous morphologies, existing detection 
methods that rely on an ellipse and a priori known radial pro¬ 
files are insufficient and new detection methods that do not have 
this limitation might provide the solution from the observational 
point of view. 

In the study of individual galaxies, the faint outer wings are 
a treasure trove for expanding our knowledge of galaxy evolu¬ 
tion. Being much less affected by internal processes, the fainter 
outskirts of galaxies preserve valuable information about the dy¬ 
namic history of the galaxy. For example, finding tidal tails 
present in these regions can be a very good indicator of possible 
previous merger events (for example, see Ellison et al. 2013). 
Accurate detection of the galaxy light profiles at large radii is 
also very dependent on the detection technique and Sky sub¬ 
traction (see Section 2.2). It has been observed that the radial 
profiles of star-forming galaxies deviate from a purely exponen¬ 
tial profile at large radii (for example, see de Grijs et al. 2001; 
Pohlen and Trujillo 2006). Successfully detecting this behavior 
in the outskirts can help explain the dependency of the star for¬ 
mation on gas surface density (Kennicutt 1989; Kregel and van 
der Kruit 2004). 

Intra-cluster light is another field in galaxy evolution that re¬ 
lies heavily on the ability to accurately detect and subtract the 
faintest parts of cluster galaxies. These regions harbor stars 
that have been stripped from the galaxies but are gravitation¬ 
ally bound to the cluster (see Presotto et al. 2014, and references 
there in). Using creative and simple hardware (Abraham and 
van Dokkum 2014), van Dokkum et al. (2015) recently found 
47 very low surface brightness but large galaxies in the comma 
cluster. Some of them could be visually confirmed in archived 
images but remained undetectable using the existing detection 
methods due to their very diffuse structure. 

Increasing the detection ability by using larger and more ac¬ 
curate hardware is one the primary reasons that new telescopes 
and cameras are being commissioned. Over the past several 
decades the instruments used for astronomical observations have 
undergone significant advancements from analog photographic 
plates to digital CCDs. However, the most commonly used detec¬ 
tion techniques still rely on the Petrosian (1976) or Kron (1980) 
approximations or de Vaucouleurs (1959) or Sersic (1963) pro¬ 
file fittings. These functions are defined on a continuous space, 
so they can be considered analog techniques. The detection tech¬ 
nique introduced in this paper exploits the digital nature of mod¬ 


ern data sets. Most of the operation is done on the noise (un¬ 
detected regions), which is also used as the basis for identifying 
false detection and segmentation results. Therefore statistically 
significant prior knowledge of signal-related properties becomes 
irrelevant. It also makes no use of any regression analysis or fit¬ 
ting which will bias the results when the targets deviate from the 
required model functional forms. 

In Section 2 the properties of the images used and some defi¬ 
nitions are provided. The new noise based solution is fully elab¬ 
orated with images showing every step in Section 3. In Section 
4 application of the concepts on a large image for increased ac¬ 
curacy and efficiency is explained. In Section 5 the success in 
detection and the purity of the proposed algorithm on mock im¬ 
ages with mock noise is analyzed. Results on mock profile and 
noise are only used for this “proof of concept” paper, a full com¬ 
pleteness, purity, number count, and etc, analysis on real images 
with no mock profiles will be provided in a companion paper 
(M. Akhlaghi et al. 2015. in preparation). Finally, in Section 6 
the usefulness of this approach for existing and future astronom¬ 
ical research is discussed. In Appendices A and B the existing 
approaches to calculating the Sky and detection and segmenta¬ 
tion/deblending are critically reviewed. Appendix C explains a 
new algorithm for finding the mode of a distribution. This algo¬ 
rithm lies at the heart of our novel sub-Sky thresholding which 
is independent of the Sky value. Appendices D and E show the 
parameter lists used for running SExtractor and NoiseChisel. 

NoiseChisel is our software implementation of the proposed 
technique in this paper. Throughout this paper, NoiseChisel is 
used to refer to the proposed new approach for detection and 
segmentation. NoiseChisel is distributed as part of the GNU As¬ 
tronomy Utilities^ which is a new collection of tools for astro¬ 
nomical data analysis and manipulation. It is being introduced 
here as the first astronomical software package that fully con¬ 
forms to the GNU coding standards^ and has been evaluated and 
endorsed by GNU. The copyright was also assigned to the Eree 
Software Eoundation to guarantee its “freedom” and facilitate 
its use and future development by the worldwide astronomical 
community similar to all their major “free” software projects. It 
integrates nicely with the GNU/Linux operating system. It can 
also be natively compiled and installed on most other major op¬ 
erating systems. The existing tools in GNU Astronomy Utilities 
were used for all the data analysis and displaying steps in this pa¬ 
per (making mock profiles, convolving, adding simulated noise, 
producing a catalog, cropping parts of large survey tiles, FITS 
image conversion to document-friendly image formats, creating 
masks, and calculating image statistics). The GNU Astronomy 
Utilities is “free” software, as defined in the GNU general public 
license (GPL),^ version 3. 

2. DEFINITIONS 

Prior to a complete explanation of the proposed new algo¬ 
rithm, some basic issues need to be addressed first. This new 
technique is non-parametric (with no functional profiles or regression- 
based fittings); therefore, demonstrating the concepts and steps 

^ https ://www. gnu. org/software/gnuastro/ 

^ The GNU coding standards define specific generally accepted requirements 
for software configuration, compilation, and installation as well as coding 
style, command line user interface, and a manual in various web-based, print, 
and command line formats. 

^ https://www.gnu.org/licenses/gpl.txt 
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Pixel intensity (c )—B{e ) 

Figure 1: Mock galaxies and their pixel distributions with and without noise. The plots have two inset images. The left insets show the image of mock galaxies after 
convolution with a PSF prior to adding noise. Logarithmic scaling has been used for (a.l) and (b.l) to show the full extent of the light prohles. The other insets all 
use a linear scale. The right insets show the image of the mock galaxies with noise added. The histogram starting from or after zero (yellow) shows the pixel count 
distribution of the non-noised mock image (left inset). The histogram is vertically truncated for a clear comparison with the noised histogram. The wider histogram 
(red) shows the histogram of the noised image (right inset). The hve vertical lines on each histogram show the median of the noised image after (4,5) a-clipping (see 
Section A.3). Insets have been truncated to be in the range of the plot’s horizontal axis: white is —250 counts and black is 1000 counts, (a) 95 random Sersic profiles 
with total magnitudes between —11 to —6 and five point sources (stars) of —12 magnitude. The zero-point magnitude for all images in this paper is set to 0.00. (b) 
A —16 magnitude Sersic profile with n = 4.00, = 30 pixels, and q = 0.25, (c) 6 large and bright Sersic prohles with —15 to —18 magnitude. The PSF is a circular 

2D Moffat prohle with FWHM= 3 pixels and p = 4.76 prior to adding noise. All the galaxies and the PSF are truncated at 5 times the effective radius and FWHM 
respectively. B is the background count, see Section 2.2. 


on actual images plays a vital role in explaining the technique 
(Anscombe 1973). In Section 2.1 a short explanation of the 
postage stamps used in this paper (including the appendices) is 
given. Fundamentally, this new technique grew out of the defi¬ 
nition of the Sky value, which is a very important concept in any 
measurement significantly affected by noise and is elaborated in 
Section 2.2. 

2. 1 . Displayed Images 

All the images in the figures are inverted: the darker the 
pixel, the larger its value. Unless otherwise stated, the displayed 
images are 200 x 200 pixels. They are all crops from original 
1000 X 1000 pixel images. For mock images, the remaining area 
is filled with pure Gaussian noise and for real images, it is from 
the actual survey. When referring to an image by histogram, the 
histogram of the number of pixels at a given count level across 
an image or section of an image is implied, similar to those of 
Figure 1. The histograms only use the pixels in the displayed 
region, not the larger image. The vertical axis values are not dis¬ 
played in all histograms beyond Figure 1, because their absolute 
values are irrelevant. 

Unless stated otherwise, the large images were used as in¬ 


put for SExtractor and NoiseChisel. The reason for setting these 
sizes was SExtractor’s modeling deficiencies in background in¬ 
terpolation and profile modeling on the corners and sides of the 
images (see Appendices A.6 and B.4). Using larger images en¬ 
sures that for both programs there is an abundant area of back¬ 
ground pixels regardless of the central object shape and size. In 
the purity analysis of Section 5.2 the pure noise regions of the 
mock images are used to identify false detections. Unless other¬ 
wise stated, all the 200 x 200 pixel cutouts are from the center 
of the larger image. 

The real images used in this paper for real demonstrations are 
either cutouts from the COSMOS survey (Scoville et al. 2007), 
Hubble Space Telescope {HST) F814W images (see Koekemoer 
et al. 2007; Massey et al. 2010), or one full CCD image taken by 
Subaru Telescope’s SuprimeCam (Miyazaki et al. 2002). 

2.2. Background, Sky, and Noise 

Throughout this paper, pixel units are assumed to be in photon- 
counts or simply ‘count’ (e~). All the methods can be readily 
applied to count sec“^ units. For example, all the displayed HST 
images were in units of counts sec“^ Let us assume that all in¬ 
strument defects—bias, dark and fiat—have been corrected and 
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the total count of a detected object, (9, is desired. The sources 
of e~ in a pixel [i^]) of the image can be written as follows: (1) 
contribution from the target object, Oij, (2) contribution from 
other detected objects, Dij. (3) undetected objects or the fainter 
undetected regions of bright objects, Vij, (4) cosmic rays, Qy, 
and (5) the background count, which is defined to be the count if 
none of the others exists on that pixel, Biy The total count in the 
pixel, Tij, can thus be written as: 

Figure 1 shows the pixel count distribution before and after 
adding noise in three sample mock images. By definition, Dtj is 
detected and it can be assumed that it is correctly subtracted, so 
that Dij can be set to zero. There are also methods to detect and 
remove cosmic rays (for example, the method of van Dokkum 
2001), enabling us to set Qj = 0. Note that in practice, Dij and 
Uij are correlated because both directly depend on the detection 
algorithm and its input parameters. Also note that no detection or 
cosmic ray removal algorithm is perfect. With these limitations 
in mind, the observed Sky value for this pixel (Sij) can be defined 
as 

Sij=Bij-hUij. ( 1 ) 

Therefore, as the detection process (algorithm and input parame¬ 
ters) becomes more accurate, or Uij 0, the Sky value will tend 
toward the background value or Sfj Bij. Thus, while Bij is an 
inherent property of the data (pixel in an image), Sij depends on 
the detection process. Over a group of pixels, for example, in an 
image or part of an image, this equation translates to the aver¬ 
age of undetected pixels. With this definition of Sky, the object 
count in the data can be calculated with 

'^ij = ^ ^ij = (^) 

Hence, the more accurately Sij is measured, the more accu¬ 
rately the count of the target object can be calculated. Similarly, 
any under- (over-)estimation in the Sky will directly translate 
into an over- (under-)estimation of the measured object’s count. 
In the fainter outskirts of an object a very small fraction of the 
photo electrons in the pixels actually belong to objects. In Fig¬ 
ure l(b.l) for example, only 19.83% of the pixels belonging to 
the central object have a value larger than 10100 counts while 
the background is 10000 ±100 counts. Therefore even a small 
overestimation of the Sky value will result in the loss of a very 
large portion of this mock galaxy. 

Based on the definition above, the Sky value is only correctly 
found when all the detected objects (Dij and Qj) have been re¬ 
moved from the data. However, as shown in Section B.l, exist¬ 
ing methods define their detection threshold using the Sky value 
(see Section B.l.2 in particular). Therefore they cannot use this 
definition in finding the Sky. They have to apply the methods 
explained in Appendix A to approximate it. The foundation of 
the technique presented here is that the detection threshold and 
initial detection process is defined independent of the Sky value 
see Figure 2 and Section 3. 

In the mock images in this paper each pixel value (Tij) is 
composed of the object contribution, Oij, and the background 5: 
Tij = B-\- Oij\ see Figure 1. The noise is added to each pixel 
by replacing its value with a random value taken from a Gaus¬ 
sian distribution with a mean of Tij and Gij = In this 



Figure 2: Flowchart of the complete detection algorithm. The steps of the top 
and bottom boxes can be seen in the example images of figures 3 and 7 re¬ 
spectively. Cl stands for convolved image, CC for connected components, .y for 
signal-to-noise ratio, A for logical ‘and’, ^ for assignment, = for comparison, 
and 0 for the empty set. Each CC is defined as a set of image pixels. I is the 
number of CCs in column 1 of Figure 7. J and K represent the number of CCs in 
the bottom and top of column 4 in Figure 7, respectively. A/, 5/, and Q are sets 
of pixels belonging to the ith connected component in their respective images. 
A, B, and C are families of sets (connected components). Di and St are the total 
signal-to-noise ratios of Bf, and Q respectively. 
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Figure 3: Demonstration of major detection steps in NoiseChisel. The rows are (a) the same mock profile as in Figure 1(b). (b) 25 mock profiles with the same 
parameters with n = 4 (see Section 3.1), (c) similar to (b), but with n = 1, (d) a real star-forming z ~ 1.5 galaxy imaged in HST/ACS F814W. The columns represent 
(1) mock profiles before the addition of noise and blank for the real image (it is blank since the signal without noise is not known in a real image; All mock images 
in this column are plotted in the logarithmic scale), (2) noisy image, (3) convolved with a 2D Gaussian kernel of FWHM= 2 pixels, (4) threshold applied on the 0.3 
quantile of the convolved image, (5) 4 connected erosions applied 2 times, (6) 8 connected opening applied 1 time, (7) estimated false detections removed (see Figure 
7; true detections are dilated two times and undetected pixels masked; pixel values are taken from the noisy image in column 2), and (8) the inverse of column 7. 
Columns 3-6 (inclusive) correspond to the top four steps of the flowchart in Figure 2. 


paper, B = 10000^“. The magnitude of each object is defined 
as —2.51og(F), were F (sum of pixel values) is in e~. In other 
words, the zeropoint magnitude is 0.00 since this paper is not 
limited to any particular instrument. So the background magni¬ 
tude is —10.00. The total magnitude is calculated from the sum 
of all the object pixels in the image (not integrated to infinity). 
All profiles are truncated at Sr^. 

3. NOISECHISEL 

The general fiowchart of the proposed technique can be seen 
in Eigure 2. In the proposed algorithm, a Sky-independent ap¬ 
proach is used to arrive at an initial detection. Individual initial 
detections are subsequently classified as true or false based on 
their signal-to-noise ratio (S/N) using the surrounding undetected 
regions as a basis. As seen in Eigure 2, this second classification 
step requires an initial measurement of the Sky value as the av¬ 
erage of the initially undetected pixels and is fully elaborated in 
Section 3.1. The substructure in each detection is then found 
and classified through defining “clumps” and “objects” which is 
explained in detail in Section 3.2. 

The input parameters'^ for NoiseChisel are specified in the 
text with — parametername=value.^ The values used for the 
parameters of NoiseChisel reported here are chosen such that 
once applied to the mock and real image tests that are shown 

^ Not to be confused with the non-parametric basis of this algorithm. By non- 
parametric it is implied that there are no functions and fittings, not that there 
are no parameters. 

^ This is similar to the long format of how the options should be called on the 
command line in NoiseChisel. 


in this paper, accurate results are achieved. The full parameter 
list can be seen in Appendix E. The criteria are to have fewer 
false detections in pure noise while detecting the faintest pixels 
of the mock profiles (see Section 5). We do not claim that these 
are “the best” set of parameters for any generic data set. 

3.1 . Detection and Sky 

The basic idea behind this detection algorithm is that if there 
is a true signal buried in the noise over a certain contiguous area, 
it will systematically (contiguously) augment that region (keep¬ 
ing the noise fiuctuations). The steps explained here are specif¬ 
ically designed to best separate these connected, augmented re¬ 
gions deep in the noise. The initial detection steps in NoiseChisel 
are schematically shown in the top group of steps of the fiowchart 
in Eigure 2 and their application to sample images is demon¬ 
strated in Eigure 3. Eirst, a very small FWHM kernel is convolved 
with the image (Section 3.1.1, column 3 in Figure 3). A very low 
threshold is then applied to the smoothed image (Section 3.1.2, 
column 4). The regions below the threshold are expanded by 
eroding those that are above it (Section 3.1.3, column 6). Since 
the objects are not completely separated yet, “opening” is ap¬ 
plied (Section 3.1.4, column 7). False detections are estimated 
and removed (Section 3.1.5) and the fainter regions of true de¬ 
tections, which were carved off during erosion are returned in a 
process of dilation (Section 3.1.6). 

Three of the examples in Figure 3 are mock and one is a 
real galaxy. Figure 3(a) is the same mock profile as in Figure 
1(b). Figures 3(b) and (c) contain 25 mock 2D Sersic profiles 
placed with equal spacing on a grid across the image. Other 
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Figure 4: Effect of convolution with wider kernels on Figure 3(c.2). The green histogram shows the count distribution of the convolved image. The blue histogram 
shows the distribution of the original image pixel count values within the green histogram’s range. The convolution kernel used is a Moffat function with p = 4.76 
(same as the PSF used for the mock images), while the FWHM is shown for each in units of pixels. Each kernel is a square of side 5 x FWHM pixels. Note how the 
convolved histogram becomes more skewed as a wider kernel is used for convolution. The range of the horizontal axis significantly decreases with increasing kernel 
width. B is the background count, see Section 2.2. 


than their magnitude and position, all of the parameters are equal 
with a position angle of 45°, an axis ratio of ^ = 0.5, and an ef¬ 
fective radius = 10 pixels. The Sersic index {n) of all the 
galaxies in each image is also the same. The bottom left pro¬ 
file is the faintest with —9.2 magnitude. The next profiles are 
successively —0.06 magnitudes brighter towards the top right. 
Such extremely faint profiles were intentionally chosen to push 
NoiseChisel (and SExtractor in Section B.1.3) to the limits. Fig¬ 
ure 3(d) shows the same steps with exactly the same parameters 
applied to a real z ^ 1.5 star-forming galaxy in the COSMOS field 
imaged with the HST Advanced Camera for Surveys (ACS) using 
the F814W filter. 

3.1.1. Convolution 

Convolution is used to maximize the ratio of an object’s peak 
signal to the local noise level (Bijaoui and Dantel 1970; Irwin 
1985). Convolution smooths the image by removing high spatial 
frequencies. This means that after convolution, “nearby” pix¬ 
els with large differences in value, which are commonly due to 
noise, will receive an intermediate value. Therefore, by smooth¬ 
ing the noise, the fraction of signal-to-noise is increased. With¬ 
out convolution, the large spatial frequencies due to noise will 
significantly limit the ability to detect a contiguous region in 
noise. Hence, in order to maximize the detection ability, con¬ 
volution is necessary. The definition of a nearby pixel and the 
functional form of how to define the intermediate value are de¬ 
fined with the kernel input into convolution. 

Figure 4 shows the noised image of 25 very faint n = I pro¬ 
files, convolved with increasingly wider (larger FWHM) kernels 
applied. As the kernel used for convolution widens, the image 
becomes more and more blurred, its histogram becomes more 
skewed, and the range of pixel values (dynamic range) in the im¬ 
age decreases with more pixels having similar values. Compar¬ 
ing the blue and green histograms in Figure 4 shows this effect. 
The skewness induced into the histogram is another manifesta¬ 
tion of the increase in the S/N of connected objects. Due to the 
skewness, separation of the brightest sections of the objects from 
the noise is facilitated. A wider convolution kernel has important 
consequences. (1) The shapes of the object tend to the shape of 
the kernel used to convolve the image. (2) The dynamic range 
significantly decreases. Therefore, separating the fainter parts of 
the brighter objects from noise becomes much harder. (3) Com¬ 
pact and faint objects, that is, those that are not wide relative to 
the convolution kernel, will be lost. 

The effects above are particularly harmful in ground-based 


images, where the FWHM is generally very large. The decrease 
in the dynamic range is not a considerable issue in the signal- 
based technique because only the brightest pixels of the objects 
are sought and the fainter parts are modeled (see Section B.4). 
As seen in the sequence of steps in Figure 3, the proposed detec¬ 
tion technique is exactly the opposite and aims to impose negligi¬ 
ble constraints and models on the detected signal. The dynamic 
range is hence extremely important to successfully find the faint 
but valuable object pixels. Therefore a sharp convolution kernel 
is used. 

Unlike contiguous (multi-pixel) signal, noise is ideally inde¬ 
pendent from pixel to pixel,^ and thus is directly connected with 
sampling or pixel size. Inverting the Nyquist sampling theorem, 
the sharpest convolution kernel can be considered an FWHM= 2 
pixel kernel and a Gaussian is used for its functional form. In all 
the examples in this paper, mock or real, ground-based or space- 
based, a Gaussian profile with FWHM= 2 pixels truncated at five 
times the FWHM is used, so that the kernel is 11 x 11 pixels^. 
The fact that the same convolution kernel gives very accurate re¬ 
sults on any type of image is one example of the objectivity of 
NoiseChisel’s output. 

Throughout NoiseChisel, the convolved image is only used 
for relative pixel values, for example in thresholding (Section 
3.1.2) or oversegmentation (Section 3.2.1). The input image is 
used for absolute pixel values, for example, the Sky value, the 
magnitude of an object, or its S/N. This distinction is not com¬ 
monly practiced in most existing techniques. For example, in 
SExtractor, the threshold value is found from the input image 
and applied to the convolved image; see Section B.1.2. 

3.1.2. Thresholding 

Pixels with values far below the inherent background or mea¬ 
sured Sky values (see Section 2.2) can harbor valuable informa¬ 
tion (see Figure 1). After accurate Sky subtraction, a real object 
will always have a positive total count. However, the standard 
deviation of the noise is always much larger than the smallest 
detectable positive count value of the real objects. For example, 
assume that the Gaussian noise standard deviation is 30^“. After 
accurate Sky subtraction, ^37% of^ the pixels that harbor 10^“ 

^ In processed images we have correlated noise which is addressed in Section 

3.1.5. 

^ The convolution kernel is created with Make Profiles which is part of the GNU 
Astronomy Utilities. The non-zero pixels in the square kernel have a circular 
shape. 

^ The cumulative distribution function of a normal distribution with mean of 
I0e~ and standard deviation of 30e~ calculated at 0e~. 
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Figure 5: Thresholds used in the demonstration examples of Figure 3 with the same labels. The blue histograms belong to the actual image (column 2 in Figure 3). 
Only the range that overlaps with the green histogram is plotted here. The green histogram is the histogram of the convolved image (column 3 of Figure 3). The thick 
curve shows the cumulative frequency of the convolved image. The thin solid vertical lines show the 0.3 quantile of the convolved image which is used to determine 
the threshold for it; this value is used to make column 4 of Figure 3. The dotted lines in the first three examples (mock galaxies) show the ideal a sky of the input image 
(blue histogram), which is I00e~. B is the background count, see Section 2.2. 


of real counts will have negative recorded pixel value. Once a 
large number of such pixels are found and the Sky value is accu¬ 
rately subtracted, the random effects of noise will be decreased 
and useful information can be extracted from such pixels. 

Therefore it is wrong to assume that the Sky value is the 
lowest possible threshold. To accurately detect the faintest pixels 
of an object, the threshold has to be less than the Sky value. This 
is very important in accurate measurement of the Sky value as 
the average of undetected pixels. 

In Noise Chisel the threshold is defined using the cumulative 
distribution of the pixels in the convolved image. To ensure that 
the threshold is below the Sky value, the quantile has to be less 
than 0.5, or the median. For the examples in this paper, the quan¬ 
tile threshold is set to — qthresh=0 . 3 , which means that 70% of 
the pixels in the convolved image will have a value above it. The 
results of applying a threshold to the convolved image can be 
seen in Figure 3, column 4, where black pixels are the ones that 
had a count above this threshold. The thresholds can be seen 
relative to the actual and convolved image histograms in Fig¬ 
ure 5. Notice that the threshold is safely below the background 
value which is 0e~ for the mock images. In Section 4.2.1 the 
complete procedure for finding the threshold on a large image is 
elaborated. 

3.1.3. Erosion 

Applying a threshold to the image (Section 3.1.2) results in a 
binary image with pixels either having a count above the thresh¬ 
old or below it. The former are known as foreground pixels 
which are black in Figure 3, column 4. The latter, called back¬ 
ground pixels, are displayed in white in the same figure. Be¬ 
cause of this extremely low threshold, it is clear from column 4 
of Figure 3 that the regions harboring true signal in them cover 
a 2D contiguous (non-porous) region of connected pixels, while 
regions where no data is present form a porous structure in the 
foreground. The figure shows that when signal is present, even 
if it is very close to the background value, it will augment the 
region that it covers. The augmentation can be seen through 
smaller and fewer white holes in the black regions that harbor 
data. The first step in Noise Chisel after applying the threshold is 
thus to exploit this porous structure and expand the holes through 
erosion of the foreground to separate the augmented regions. 

Erosion is an operation in morphological image processing 
(Starck and Murtagh 2006; Gonzalez and Woods 2008). By 
eroding the foreground, it is implied that all foreground pixels 


that neighbor a background pixel are carved off^ the foreground 
and changed to a background pixel. The neighborhood of, or 
connected pixels to, a pixel is defined based on a structuring el¬ 
ement and generally has two types: 4 and 8 connectivity (see 
Figures 6(a) and (b)). Eroding based on 8 connectivity will re¬ 
move more pixels. Therefore, in this step, NoiseChisel uses a 4 
connected structuring element to erode the foreground. The sta¬ 
tistical properties of this and other operations in mathematical 
morphology are reviewed in Dougherty (1992). 

Eroding the foreground expands the holes (white regions in 
Eigure 3 column 4) and connects them to each other. The best 
number of erosions (— erode) and the type of connectivity (— erodengb) 
applied to an image is intertwined with — qthresh. The results 
of respectively setting them to 2 and 4 can be seen in Eigure 3, 
column 5. If a set of foreground pixels had the shape of the non¬ 
white pixels of Eigure 6(c), after applying this erosion the pixels 
with the two brighter shades of gray would be carved off and 
only the central 9 pixels would remain. 

3.1.4. Opening 

Dilation is the inverse of erosion in morphological image 
processing. When dilating the foreground, if a background pixel 
touches a foreground pixel, it will be changed to the foreground. 
Applying these two steps in order, namely, first eroding then di¬ 
lating the foreground, is known as opening. 

Opening is very useful for separating regions of the fore¬ 
ground that are only connected through a thin thread of pixels. 

It also has the benefit of removing extremely small objects (see 
Vollmer et al. 2013, for such an application in astronomy). If the 
image is first eroded n times and then dilated the same n times, 
then n is referred to as the “depth” of the opening. Note that 
the erosion discussed here is part of opening and separate from 
the one in Section 3.1.3. The depth can be set in NoiseChisel 
through the parameter —opening and the type of connectivity 
can be set through — openingngb. In the examples shown in Eig¬ 
ure 3, they are set to l and 8, respectively. The opened image 
can be seen in column 6 of Eigure 3. Compared with column 5, 
it is clear that, as expected, opening has separated the foreground 
into separately connected regions or islands without bridges. 

Multiple erosions and dilations with the same structuring el¬ 
ement will result in the outer shapes of the foreground becoming 
like a diamond in 4 connectivity or a rectangle in 8 connectiv¬ 
ity. Therefore, since the erosion of Section 3.1.3 was based on 

^ Hence the name NoiseChisel: a tool for carving off noise, similar to what 
wood chisels or stone chisels do on their respective materials. 
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a b c 

Figure 6: Structuring elements. The black pixel is the pixel under consideration 
and the gray pixels are the neighbors of that pixel, (a) 4 connected neighbors, 
(b) 8 connected neighbors, (c) Non-white pixels show the smallest area an object 
should have to be detected with the given two 4 connected erosions (see Section 
3.1.3) and one 8 connected opening (see Section 3.1.4). It is assumed that all the 
pixels outside the box shown are background (white). The three shades of gray 
show the pixels that are removed in each step (brightest first). After the three 
steps finish, the inner 9 pixels remain. In processing, white pixels have a value 
of 0 and others have a value of 1. 

a 4 connected structuring element, in this step an 8 connected 
structuring element was used to compensate for the shape loss 
in the previous step as well as more successfully removing thin 
connections. The composite effect of the two steps of erosion 
in Section 3.1.3 and the opening done here can be seen in Fig¬ 
ure 6(c). Before the dilating section of the opening step, only 
the central black pixel will remain. If a region above the thresh¬ 
old is any smaller than the colored pixels of Figure 6(c), it will 
disappear after this step since no pixels remain for dilation. 

The erosion and dilation defined here are based on a binary 
image as discussed in Section 3.1.3. It is also possible to de¬ 
fine such operations on grayscale images where a pixel can have 
a range of values, not just 0 or 1. Ferret et al. (2009) used a 
combination of such grayscale erosion and dilation (hit-or-miss 
transform) for a demonstration on detecting low surface bright¬ 
ness galaxies. Since the structuring element is no longer a bi¬ 
nary, the variation of pixel values over the structuring element 
becomes very important and has to be modeled based on spe¬ 
cific elliptical parameters, surface brightness profiles, and the 
PSF; see Figures 4 and 5 of Ferret et al. (2009). Therefore, it is 
optimal only for profiles that satisfy the count distribution of the 
structuring element. Applying it generically to any image with 
any distribution of galaxy morphologies will require a high level 
of customization that may fail for complex galaxy morphologies 
like the real galaxies shown in this paper. 

3.1.5. Defining and Removing False Detections 

The foreground is now composed of separately connected 
objects that we define as initial detections. However, on the 
mock images (particularly comparing Figure 3(a.l) with (a.6)), 
it is clear that some detections can be associated with true signal 
while most cannot. These “false” detections are sets of pixels 
with random values (above the very low threshold on the con¬ 
volved image) that were randomly positioned close enough to 
each other to pass the erosion and opening steps of detection 
explained above. As the threshold decreases, the probability of 
such random positioning increases. If a larger threshold, number 
of erosions, or number of openings is used, false detections will 
decrease or completely disappear. 

The problem with more stringent detection parameters is that 


besides rejecting more inherently faint objects, they will also re¬ 
move the valuable fainter parts of the final detections. Therefore, 
a classification scheme is defined to find and remove false detec¬ 
tions very accurately. This allows for keeping the valuable faint 
pixels of bright detections and most of the detections that would 
have been lost otherwise. 

The bottom group of steps in the flowchart in Figure 2 show 
all the steps used to define and remove false detections. The 
same steps can also be seen in practice in Figure 7 and elaborated 
here. The detections have divided the pixels into 2 sets: (1) The 
undetected sky, Rs (white regions of column 1 in Figure 7) and 
(2) the detected regions, Rd (black regions in column 1 of Figure 
7). It can be reasonably assumed that no significant contribution 
from a detectable object exists in Rs. Some of the detections in 
Rd harbor signal and some are localized random noise pixels, 
namely false detections. 

In order to identify the false detections m Rd, a secondary 
detection process is applied to the pixels of both Rs and Rd inde¬ 
pendently. The detections in the former will provide a scale on 
which the certainty of the detections in the latter can be defined. 
Let Sa and be the average and standard deviation of the count 
in Rs, respectively. The second threshold is defined by these two 
parameters. Note that this Sa is not the final Sky value because 
some of the high valued, localized noise pixels (false detections) 
have systematically been removed from it; it is just used here as 
an approximation. 

The user is free to set this second threshold based on SaF 
dthresh XG^. In the examples here, — dthresh=-0.1. The re¬ 
sult of applying this threshold to the detected and blank regions 
can be seen in Figure 7, column 2. Note that as discussed in 
Section 3.1.1, absolute values such as Sa and G^ come from the 
input image. Therefore, while the pixel indices comprising Rs 
and Rd are defined from column 6 of Figure 3, which used the 
convolved image, the pixel values for finding and applying this 
second absolute threshold come from the actual image. The pix¬ 
els that were above this threshold are marked in black in column 
2 of Figure 7. 

As discussed in Section 3.1, a significant advantage of data 
over noise detections is that the data augment a contiguous re¬ 
gion. In order to exploit this contiguous property of true detec¬ 
tions, all the holes that are engulfed in 4 connected foreground 
regions will be filled. Note that holes are white pixels in Figure 
7, column 2. The result of filling such holes can be seen in col¬ 
umn 3 of Figure 7. Such holes are most probably due to noise 
since the 4 connectivity of all the surrounding pixels is a very 
strong constraint (compare Figure 7(b.2.2) and (b.3.2)). Com¬ 
bined with the next step (opening), filling holes in this manner 
significantly enhances the ability to identify a faint contiguous 
signal. 

To remove the thin connections between thicker regions and 
obtain individual connected components, one level of 4 con¬ 
nected opening is applied. The result can be seen in column 4 of 
Figure 7. The two images of column 4 in Figure 7 are now com¬ 
posed of separate connected components (pseudo-detections) that 
were created in exactly the same manner but on different regions 
of the input image. The parameter used to quantify the defini¬ 
tion of a false detection is the total S/N (S/Nr) of each pseudo¬ 
detection which is a combination of its total count and area. 

Let F be the average count in each pseudo-detection and 
N be its area. The areas vary widely, ranging from one pixel 
to thousands of pixels. Therefore in order to get a reasonable 
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Figure 7: Identifying and removing false detections, demonstrated for three of the examples from Figures 3(a), (b) and (d). Column 1 is the same as column 6 in Figure 
3. Columns 2 to 4 apply the same steps (explained below) on the undetected regions (top) and detections (bottom) separately. Column 2: Applying the threshold. 
Column 3: Fill holes surrounded by one 4 connected region above the threshold. Column 4: Application of opening to the previous column and small regions removed. 
Column 5, top: distribution of the total signal-to-noise ratio of all the dark connected components in column 4 (top). The dashed vertical line in the histogram shows 
the —det quant=0.99 quantile of this distribution. This value is used as a threshold to define false detections in column 4 (bottom), (a.5.1) is from a nearby large mesh 
(see Section 4) because the number of connected components in (a.4.1) is less than —minnumfalse=100. Column 5 (bottom): all detections with total signal-to-noise 
ratio above that threshold. The arrows in (a) visualize the steps which are schematically shown in the lower box of Figure 2. 
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average count comparison, all pseudo-detections with an area 
smaller than — detsnminarea are removed from the analysis in 
both Rs and R^. In the examples in this section, it is set to 
15. Note that since a threshold approximately equal to the Sky 
value is used, this is a very weak constraint. For each pseudo¬ 
detection, S/Nr can be written as, 

vny-s.) ^ 3 ^ 

^NF+Ndl ^F + al 

See Section 3.3 for the modifications required when the input im¬ 
age is not in units of counts or has already been Sky subtracted. 
The distribution of S/Nr from the objects in Rs for the three ex¬ 
amples in Figure 7 can be seen in column 5 (top) of that figure. 
Image processing effects, mainly due to shifting, rotating, and 
re-sampling the images for co-adding, on the real data further 
increase the size and count, and hence, the S/N of false detec¬ 
tions in real, reduced/co-added images. A comparison of scales 
on the S/N histograms between the mock ((a.5.1) and (b.5.1)) and 
real (c.5.1) examples in Figure 7 shows the effect quantitatively. 
In the histograms of Figure 7, the bin with the largest number 
of false pseudo-detections respectively has an S/N of 1.89, 2.37, 
and 4.77. 

The S/Nr distribution of detections in Rs provides a very ro¬ 
bust and objective scale that can be used to select or reject a given 
detection in Rd. The level of certainty for true detections can be 
input by the user and is defined based on — detquant which 
is the desired quantile in the S/Nr distribution of Rs. In the 
flowchart shown in Figure 2, this is displayed with Dj > Sq. The 
accepted detections from the second thresholding can be seen in 
the bottom row of column 5 in Figure 7. Any initial detection 
that overlaps with at least one of the accepted pseudo-detections 
{Bj n A; 9 ^ 0 in Figure 2) is defined as a true detection, while 
those that don’t are considered a false detection and removed. 
The successful detections can be seen in the last two columns of 
Figure 3 (after a final dilation; see Section 3.1.6). 

Outside of astronomical data analysis, the technique pro¬ 
posed here to separate true from false detections (and later clumps; 
see Section 3.2.1) can be considered a form of anomaly detec¬ 
tion in data mining. An anomaly or outlier is defined as “an ob¬ 
servation which deviates so much from other observations as to 
arouse suspicion that it was generated by a different mechanism” 
(Hawkins 1980, page 1). This is a very general and qualitative 
definition applicable to any form of anomaly (see Chandola et al. 
2009, for a review). 

The S/N of the connected components in Rs (after applying 
the threshold, filling holes, and opening) gives a scale to define 
the outliers in Rd (true pseudo-detections). However, a quantile 
and not the histogram is used (the histogram is only displayed 
in Figure 7 for demonstration). Therefore based on the classifi¬ 
cation of Chandola et al. 2009, the true/false classification tech¬ 
nique proposed here can be considered a non-parametric, “quan¬ 
tile based” statistical anomaly detection technique. 

This technique for determining an “anomaly” (or a real ob¬ 
ject within noise) has some similarities to the Higher Criticism 
statistic (Donoho and Jin 2004). It was used in detecting non- 
Gaussianities in the cosmic microwave background (Cayon et 
al. 2005) with the modification of using p-values. However, the 
higher criticism statistic is derived from and applied to the same 
set of observations, while here the S/N threshold is found using a 
completely separate set of pixels (those in Rs). It is then applied 


to another set of pixels {Rd). Also, while the Gaussian/Poisson 
distribution is assumed for the noise, there is no such assumption 
for the distribution of the S/N. 

3.1.6. Final Dilation 

The final step is to restore, through dilation, the pixel layers 
of each remaining initial detection that were removed by erosion 
in Section 3.1.3. Dilating — erosion times would presumably re¬ 
store the object pixels out to the initial — qthresh quantile. The 
user can set the number of final dilations with the —dilate op¬ 
tion. Images of astronomical objects never have a strong cutoff. 
Even if they have a physical cutoff, the PSF will soften any sharp 
boundary except cosmic rays which are independent of the PSF. 
Therefore it can be assumed that the object extends beyond the 
initial threshold. To exploit this fact, this final dilation is only 
based on 8 connectivity and in the examples of this paper; it is 
set to — dilate=3. 

3.2. Segmentation 

The detection thresholds used are extremely low in NoiseChisel 
(see Section 3.1.2). Combined with the fact that galaxies have 
no clear cutoff, it often happens that several apparently nearby 
objects on the image will be detected as one region. There are 
two approaches to separating a detected region into potentially 
several sub-detections or to find substructure: segmentation and 
deblending. 

Segmentation is the procedure that assigns each pixel to one 
sub-component. Deblending, on the other hand, identifies the 
contribution in each pixel from different blended sources. There¬ 
fore in deblending, each pixel can belong to more than one ob¬ 
ject. Deblending is more realistic: the sum of the light profiles of 
multiple overlapping objects specifies the final pixel value. How¬ 
ever, to be accurately done, deblending will require parametric 
analysis and fitting algorithms. Segmentation, on the other hand, 
is a very low-level measurement and only goes so far as to as¬ 
sign each pixel to the sub-component that contributes most to it. 
Therefore segmentation provides the best starting point for the 
higher level deblending if it is needed. NoiseChisel will not do 
deblending but only performs segmentation. 

In order to segment each detection into sub-components, clumps 
(Section 3.2.1) and objects (Section 3.2.2) are defined. A de¬ 
tected region might contain one or more objects and anything 
from none to multiple “true” clumps. Finding the clumps in an 
object or detected region is a low-level measurement that pro¬ 
vides us with a multitude of false clumps mixed with possibly 
true clumps. True clumps can be found robustly by studying the 
noise properties in the vicinity of the detected regions. This pro¬ 
cess is very similar to how true/false detections were separated 
in Section 3.1.5. In Section 3.2.1 the proposed definition for a 
clump and the true/false classification is explained. The segmen¬ 
tation of a detected region to objects is a high-level measurement 
thoroughly discussed in Section 3.2.2. 

3.2.1. Finding “True” Clumps 

Pixels at a local maximum are defined as those whose value 
is larger than their 8 connected neighbors. The 8 connected re- 

We will be working on such a fitting tool for deblending as part of the GNU 
Astronomy Utilities to start using the outputs of NoiseChisel. 
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Figure 8: Pixel-by-pixel demonstration of oversegmentation (see Section 3.2.1 for a complete explanation of the algorithm). The order is from the top left to the 
bottom right. First: a small region in the center of Figure 9(b.2). The rest show how pixels are labeled in order of their value. The shades of gray become darker as the 
labels increase. The darkest shade of gray represents river pixels. 


gion around each local maximum whose pixels all have a lower 
count than the peak is defined as a clump. Note that the clump 
can be much larger than 9 pixels. 

The approach introduced here for finding the clumps associ¬ 
ated with all the local maxima of an image is conceptually very 
similar to and inspired by the method proposed by Vincent and 
Soille (1991), without using layers because of the biases they 
produce (see the discussion in Section B.2). Figure 8 shows the 
process for 124 pixels of the central region of Figure 9(b.2). 

A zero valued array with the size of the image is first created 
to store the labels of the different clumps. All the pixels of each 
object are ordered based on their value. Starting from the bright¬ 
est pixel and in decreasing order, the pixels are labeled based on 
the following conditions. 

1. If the pixel is not 8 connected to any other labeled pixels, it 
is a local maximum and is assigned a new label. This can 
be seen as new labels (darker shades of gray) are added in 
Figure 8. 

2. If its 8 connected neighbors all have the same label, the 
pixel is given that label as well. In Figure 8, this occurs 
most often (when the colored areas expand). 

3. If its neighbors have different labels, like a river flowing 
between two mountains, the pixel is a local minimum. 
The river pixels act as separators of each mountain or 
clump. River pixels are the darkest pixels in Figure 8. The 
pixel-by-pixel construction of the rivers can be seen from 
the middle of the third row in Figure 8 as the growing 
regions start touching. In this setup the rivers are 4 con¬ 
nected. If 4 connectivity was chosen for checking neigh¬ 
bors, the rivers would be 8 connected. 

The application of this simple algorithm to the non-detected 
regions of the image can be seen in column 3 of Figure 9 and 

As Vincent and Soille (1991, page 583) put it, “In the field of image process¬ 
ing and more particularly in mathematical morphology, grayscale pictures are 
often considered as topographic reliefs. In the topographic representation of a 
given image /, the numerical value (that is, the gray tone) of each pixel stands 
for the elevation at this point. Such a representation is extremely useful since 
it allows one to better appreciate the effect of a given transformation on the 
image under study.” In this analogy, the noisy image can be considered as a 
mountain range and when rain comes, the water gathers in the local minima 
of those mountain ranges to form rivers. 


on the detected regions in column 5. Similar to our discussion 
in Section 3.1.5, the region of no detection is used as a basis 
or calibrator to find the true clumps over the detected regions. 
This algorithm for finding the clumps uses their relative values. 
Therefore as discussed in Section 3.1.1, the convolved image 
is used to find the index of pixels in each clump. The small 
convolution kernel is beneficial for this step because the spatial 
resolution is very important for an accurate result. This pixel- 
by-pixel expansion process provides the ultimate usage of the 
dynamic range of the convolved image. 

A very large fraction of the clumps found over a detected ob¬ 
ject (column 5 in Figure 9) are due to background and correlated 
noise. To identify the false clumps, the segmentation results on 
the noise regions of the image (column 3 of Figure 9) are used as 
a reference. As discussed in Section 3.1, it can be assumed with 
high certainty that there is no significant signal (from a true phys¬ 
ical object) in the noise regions. Furthermore, the noise regions 
in the vicinity of the object would contain the same undetected 
astronomical objects and instrumental and data analysis biases 
as the object pixels. Therefore this region is the most accurate 
and objective reference available to judge between a true and a 
noise clump. 

On the original unconvolved and not Sky subtracted image, 
Fi and Fq are respectively defined as the average count in a clump 
and on all river pixels surrounding it. The average value of the 
river pixels around each clump approximately show the base el¬ 
evation or count if the clump was not present. Let Ni be the total 
area or number of pixels inside a clump. The S/N of each clump 
can be written as, 

^/w^+^K ^/fr+^v 

See Section 3.3 for how this equation can be modified for 
images that are not in units of counts or have already been Sky 
subtracted. The S/N distribution of all the clumps in the blank 
Sky can be seen for all the examples of Figure 9 in column 4 of 
that figure. The — segquant=0.99 quantile of the distribution is 
used to find the S/N threshold to accept or reject a clump over the 
detected object. This threshold is thus completely independent 
of the particular object in which the potential clumps are embed¬ 
ded. Note that the criteria in SExtractor and some other existing 
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Figure 9: Finding clumps through oversegmentation. Column 1: input images. Column 2: convolved image (Section 3.1.1). Column 3: oversegmentation applied 
to non-detected regions on the convolved image. Column 4: signal-to-noise ratio histogram (Equation (4)) of all segments in column 3 with an area larger than 
—segsnminarea=25. The dashed line shows the position of the — segquant=0.99 quantile. Column 5: oversegmentation applied to the central object. Column 6: 
clumps in all the detected objects with signal-to-noise ratio larger than the threshold of column 4. Note that some detections have no true clumps. 


tools for deblending an object is the fraction of counts in each 
clump to that in the parent detection (see Section B.2 for a dis¬ 
cussion of the resulting biases). Since no layering is necessary 
and the clump counts are not compared to the object anymore, 
the clumps in various galaxies can now be objectively compared 
with each other. 

Since the average count in each clump is very important 
in the S/N calculation, only clumps having an area larger than 
—segsnminarea are considered in the comparison. Any clump 
that is smaller than this area will be discarded as noise. For 
the examples in this paper it is set to 25. This is larger than 
—detsnminarea, because the areas there were found by apply¬ 
ing a threshold to the input image, while the clumps here are 
found on the convolved (smoothed) image where the areas be¬ 
come larger. Clumps smaller than this area add very strong 
skewness to the S/N distributions. 


3.2.2. Object Segmentation 

The very low thresholds that were used (see Section 3.1.2) 
cause objects that are close enough to each other on the image to 
be blended in one detection. Three real examples are displayed 
in the first column of Figure 10. The basis for separating po¬ 
tential objects in a given region is the true clumps which have 
been found in that detection. Therefore, if a detected region has 


none or only one clump, it is considered to be only one object. 
If there is more than one clump in a detected region, the clumps 
are grown until a certain threshold of the input image. If the 
Sky value (average of non-detected pixels) and its standard devi¬ 
ation are labeled with S and a^, then the threshold to stop clump 
growth is 5'+gthreshxa5. In the examples of Figure 10 it is set 
to — gthresh=0.5. 

The process of growing the clumps is very similar to the al¬ 
gorithm in Section 3.2.1 (Figure 8) with the exception that no 
new labels are added. If a pixel has no labeled neighbors, it is 
kept in a queue to be checked on a next loop. The loop continues 
until no new pixel can be labeled. The grown segments of each 
clump of the examples in Figure 10 can be seen in column 3. 

Objects are defined based on the average S/N of the river 
pixels between the grown clumps of Figure 10, column 3. Two 
grown clumps are defined as separate objects if the river between 
them is below a user-defined S/N. If it is larger, then the connec¬ 
tion between the grown clumps is too strong to be regarded as 
separate objects. Therefore the grown clumps are considered as 
parts of one object. Let Fij represent the average count on the 
river between the grown clumps i and j. The average S/N of the 
river between these two clumps is defined as: 


S/No- 


Fjj-S 
a / ^ij + 


( 5 ) 
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Figure 10: Segmenting objects. Column 1: actual images. Column 2: true clumps in the detections (see Section 3.2.1). Each clump is color-coded based on its label. 
The underlying light gray region shows the detection region. Column 3: the clumps are grown until — gthresh=0 .5 (see Section 3.2.2). Column 4: all objects with 
boundary signal-to-noise ratio larger than — ob jborder sn=l are given the same label. Column 5: the remaining pixels of the region are filled by growing the labels of 
column 4. In the last column only the central segmented objects are shown for clarity. 


If S/Nij <ob jbordersn, then these two grown clumps are con¬ 
sidered separate objects and if not, they are defined as belong¬ 
ing to the same object. In the examples of Figure 10 this is 
set to — objbordersn=l. S/Nij{=S/Nji) is calculated for all the 
grown clumps of each detection. Finally, all the clumps that are 
connected to each other within a detection are given one label 
(see column 4 of Figure 10). Note that if two grown clumps 
have S/Nij <ob jbordersn, and they are both connected to a third 
clump with S/Nik >objbordersn and S/Njk >objbordersn then 
all three grown clumps will be given the same label. The rivers 
between two grown clumps are only one pixel wide. Therefore 
in the calculation of Ftj, the count of each pixel is taken as the 
average of that pixel and its 8 connected neighbors. Hence, in 
practice the rivers between two grown clumps are 3 pixels wide. 

Because of the absolute nature of the two parameters for 
object definition (as opposed to the relative nature in the case 
of detection and clumps), namely the user directly providing 
—gthresh and —objbordersn, the object definition is not ro¬ 
bust and objective like that of detections and clumps. In fact 
—ob jbordersn is the only user given S/N value in NoiseChisel. 

The definition of objects made here is purely based on the 
one image input to NoiseChisel. The example above shows that 
this definition of objects is fundamentally subjective (directly de¬ 
termined by the particular data-set and the user and not neces¬ 
sarily based on physical reality). Two galaxies at vastly different 
redshifts might be sufficiently close to each other, on the same 
line of sight, that they are each defined as clumps in one larger 
object. However, due to surface brightness limits, the connection 


(for example, through a spiral arm or a tidal tail) of a real star¬ 
forming region of a galaxy might not be detected. Therefore that 
region will be identified as a separate object. This also applies 
to rings and filaments. In order to solve these issues, the infor¬ 
mation in one data set (image in this case) is not sufficient, and 
ancillary data, for example, from spectroscopic data or images 
in other wavelengths, will be necessary. Therefore the physical 
reality of the ‘objects’ defined here is beyond the scope of any 
analysis based solely on one image. 

3.3. S/N Equation Modifications 

The standard deviation of the non-detected regions (a^), in¬ 
corporates all sources of error: Poisson noise, read-out noise, 
etc. Hence, if the image is in units of counts and the backgrounds 
are not already subtracted once, then Equations (3) - (5) can be 
used no matter if the noise is background-dominated or read- 
out-noise-dominated. However, input images do not necessarily 
satisfy this condition. For example, the HSTf ACS images shown 
above are processed (already Sky-subtracted) and distributed in 
units of counts per second. 

When the image is in units of counts s“^ and if the noise is 
background-dominated, all the terms in the equations have to be 
multiplied by a constant in units of time. In finding true detec¬ 
tions and clumps that depend on the quantile of the S/N distribu¬ 
tion, the absolute value of this constant makes absolutely no dif¬ 
ference. However, in identifying objects, it is important, because 
the user specifies an absolute S/N value for rivers. NoiseChisel 
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will automatically detect if an image is in units of counts or 
counts s“^ using the minimum standard deviation in all the meshes 
(see Section 4). If this minimum is larger than unity, then it is 
assumed that the image is in units of counts. If not, the inverse 
of this minimum (which has units of time) will be used as the 
time constant discussed above. 

If the input images are processed, then Sky subtraction has 
already been applied. Therefore the error in the average count F 
in the denominators should be F + aj, not F alone as is currently 
present in the equations. It is up to the user to specify this condi¬ 
tion through the — skysubtracted parameter (see the manual for 
more details). 

4. LARGE (REAL) IMAGES 

Real astronomical images are not as small and clean as the 
examples in the figures in Section 3. They are often much larger 
or might contain masked pixels. The former often causes varia¬ 
tions in background and noise (see Eigure 11). These gradients 
do not necessarily have to be created by the background, but they 
might be due to processing, for example, bad fittings in fiat field¬ 
ing. Even if the image is very clear, and the targets are already 
accurately defined, the statistics when using such small postage 
stamps will not be accurate. By statistics we mean the number 
of false detections and clumps in the non-detected region that 
go into calculating the S/N quantile; see Section 3.1.5 and 3.2.1. 
Therefore much larger postage stamps, for example, from sur¬ 
veys, have to be used. As discussed in Section 2.1, the input 
FITS images of all the figures so far were much larger than those 
shown here. 

The image is considered to be a grid of meshes. All the steps 
explained in Section 3 can be applied independently on each 
mesh. Since the operation on each mesh is completely inde¬ 
pendent of the rest, parallel processing can be applied to signif¬ 
icantly improve processing time. This is only possible because 
all the methods in Noise Chisel are non-parametric. 

Some procedures are directly linked to the potential gradi¬ 
ents in the image, for example, convolution (Section 4. 1), thresh¬ 
olds (Section 4.2.1), and Sky subtraction (Section 4.2.2), while 
others are only defined after the gradients have been removed 
and need a large area for statistical accuracy, for example, re¬ 
moving false detections (Section 3.1.5) and finding clumps (Sec¬ 
tion 3.2.1). To be able to adequately do both jobs, two mesh sizes 
are defined: a small mesh specified by — smesh for the former 
class of operations and a large mesh specified by —lmesh for the 
latter. Each specifies the sides of a square mesh. As explained in 
Section 2.1, in the examples of Section 3, — lmeshsize=200 (so 
one mesh covers the displayed postage stamps) and for Eigure 11 
it is set to 256. In all the examples of this paper, — smeshsize=32. 

4.1. Convolution 

Image convolution is mostly done using the discrete Eourier 
transform in the frequency domain (see Gonzalez and Woods 
2008). When the convolution kernels are large, for example, 
when the image PSF is used, this technique provides signifi¬ 
cant performance benefits. However, it fails on the edges of the 
image. Spatial domain convolution, on the other hand, can be 
corrected to account for the edge effects or masked pixels and 
when the kernels are small it can be even faster. The kernels in 


NoiseChisel are small so added to its extra capability, spatial do¬ 
main convolution is used (see the Convolve section of the GNU 
Astronomy Utilities manual for a complete explanation). 

4.2. Interpolation and Smoothing 

Some meshes will not be able to provide an input in the 
grid, for example, because of a large object that is larger than 
the mesh. Therefore interpolated values will be used for unsuc¬ 
cessful meshes, and finally, the mesh grid is smoothed. We will 
discuss the case for the quantile threshold (Section 4.2.1) and 
Sky value (Section 4.2.2) more accurately here. 

Following the non-parametric nature of NoiseChisel, func¬ 
tional interpolation techniques like spline or cubic interpolation 
will not be used. They can result in strong outliers on the edges 
or corners like those in Figure 16(a)-(c) and (1). To interpolate 
over each blank mesh, the median value of the nearest — numnearest 
acceptable meshes is used. The nearest non-blank neighbors are 
found efficiently through a breadth-first search strategy in graph 
theory. In such a median interpolated grid, the result will not be 
smooth. Hence an average filteris applied to the interpolated 
grid to make smoother variations between neighboring meshes. 
The final interpolated grid for the Sky value can be seen in Fig¬ 
ure 11(c). 

4.2.1. Threshold 

The quantile threshold of Section 3.1.2 can only have a com¬ 
parable value between all meshes, if the signal in the mesh is not 
significant compared to the noise. Otherwise, with more signal 
contributing to the mesh, the pixel distribution in the mesh will 
become more skewed and thus the quantile value found will also 
shift and not be comparable between different meshes. In order 
to find the threshold, the novel algorithm of Appendix C is used 
to find the mode of the image. The mode acts as a gauge for the 
amount of data that is mixed with the noise in that mesh. 

With more signal contributing to the mesh, the average, me¬ 
dian, and mode will shift (with decreasing rates) to the positive 
(see Figure 1 and Appendix C). Exploiting this property, the sig¬ 
nificance of signal in noise of a mesh regardless of its morphol¬ 
ogy can be assessed using the difference between the median and 
the mode such that the closer the mode is to the median, the less 
significant signal there is in the mesh. Through the parameter 
—minmodeq, the user can set the minimum acceptable quantile 
for the mode. In this paper it is 0.49. Note that the skewness 
induced by convolution, further skews the distribution, or, dis¬ 
tances the mode and median (see Section 3.1.1 and Figure 4). 
Hence the 0.4 9 value is a very strong constraint. 

4.2.2. Sky and a sky 

As defined in Section 2.2, the Sky value is the average of 
undetected pixels in each mesh and its error is their standard 
deviation. Figure 11(a) shows an example Subaru Telescope 
SuprimeCam image that has undergone initial processing after 
fiat fielding and prior to Sky subtraction. The visible gradients 
in the input image are most probably due to problems in fiat 
fielding rather than actual variation in the Sky. Finding the cause 

The value of each mesh is replaced with the average of it and its 8 connected 
neighbors (see Figure 6(b)). 
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Figure 11: Applying NoiseChisel on a raw image prior to Sky subtraction from Subaru SuprimeCam. The image is 2048 x 4177 pixels, (a) The raw image after bias 
subtraction, dark subtraction and flat fielding have been applied, (b) All the detected regions in the image have been removed to expose the Sky. (c) The interpolated 
and smoothed Sky value defined as the average of undetected regions over a mesh grid, (d) The Sky value on each mesh is subtracted from the respective region of the 
image. Notice how the gradients in the input image have been mostly removed, (e) Enlarged image of white box in (a), (f) Zoom-in of the same region in (b). The 
color truncation values in (a)-(c), (e), and (f) are exactly the same; for (d) they were set such that the bright star has the same apparent diameter. The white pixels in (a) 
and (d) are masked pixels. 


of this gradient is very important, but beyond the scope of this 
paper. No matter how it was created, a good detection algorithm 
has to be as resilient to such gradients as possible. This particu¬ 
lar image was intentionally chosen to demonstrate NoiseChisel ’s 
ability to account for such gradients. 

Figure 11(b) shows the detected objects removed from the 
input image. The image is covered with a mesh grid of size 
—smesh. The Sky on each mesh is found from the average of the 
undetected pixels in that mesh. Only those meshes with a suffi¬ 
ciently large fraction of undetected pixels are used. If too much 
of the mesh area is covered with detected objects, that mesh will 
not be used because of the potential fainter wings penetrating to 
the undetected regions; compare Figures 3(a.l) and (a.8). The 
fraction can be set by — minbf rac, in this paper it is set to 0 .5. 

Cosmic rays can be a problem for determining the average 
and standard deviation. For example, long exposures in raw 
HST/ACS images are heavily peppered with cosmic rays (see Fig¬ 
ure 15 for a normal example). Since cosmic rays are very sharp 
and can be very small, some will not be detected by the detec¬ 
tion algorithm. Therefore a simple calculation of the average and 
standard deviation will be significantly biased if cosmic rays are 
present. To remove the effect of cosmic rays, a-clipping with 


the same convergence-based approach as S Extractor is used (see 
Section A. 3). Note that by this step the mode and median are 
approximately equal, so that a-clipping is just used to remove 
the effect of cosmic rays on the average and standard deviation. 

The interpolated and smoothed (see Section 4.2) Sky map 
can be seen in Figure 11(c), each pixel in this image represents 
a 32 X 32 pixel mesh. Finally in Figure 11(d), the Sky value 
for each mesh is subtracted from the pixels in that mesh. The 
gradients and differences between the different CCD chips have 
been significantly suppressed. Notice that even with this small 
mesh size some of the strong top right gradient still remains. 

5. ANALYSIS 

NoiseChisel is based on detecting an object deeply buried in 
noise using its fainter parts (see Section 3). In contrast, existing 
techniques rely on detecting the brightest regions and growing 
those (based on a priori models) to include the fainter parts (see 
Appendix B). Hence these methods are fundamentally different. 
The detection accuracy is analyzed here using three fundamental 
and basic statistics that can be done with mock profiles in mock 
noise. Comparing the Sky and background value when the latter 
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Figure Label 

Median 

Mode 

a-clip 

Converge 

a-clip 

(4,5) 

SExtractor 

Average of 
Undetected 

NoiseChisel 

Average of 
Undetected 

NoiseChisel 
Nf (Section 5.2 

12(a) and 18(a) 

5.6 

5^ 

5.4±99.9 

5.6±99.9 

5.1±99.9 

2.2±99.2 

3 

12(b) and 18(b) 

6.6 

6.6 

6.4 ±100.7 

6.6 ±100.7 

5.9 ±100.7 

3.5 ±100.1 

4 

12(c) and 18(c) 

4.8 

2.2 

5.6± 101.2 

4.8± 101.2 

3.9 ±100.5 

3.1 ±100.4 

2 

12(d) and 18(d) 

5.9 

5.9 

6.3 ±102.2 

5.9 ±102.2 

3.1 ±100.4 

2.2± 100.1 

2 

13(a) and 19(a) 

3.4 

1.3 

5.1 ±29.0 

3.1 ±28.6 

-0.9 ±25.0 

-1.2±25.0 

N/A 

13(b) and 19(b) 

2.0 

0.7 

3.1 ±27.9 

1.9±27.8 

-0.7 ±25.8 

-0.8 ±25.9 

N/A 

13(c) and 19(c) 

2.1 

0.7^^ 

2.5 ±22.0 

1.6±21.7 

-0.9±19.7 

-1.1±19.8 

N/A 

13(d) and 19(d) 

4.1 

3.0 

5.2±28.6 

4.1 ±28.5 

-1.1±25.5 

-1.2±25.6 

N/A 

13(e) and 19(e) 

38.9 

32.3 

41.1±115.8 

36.1 ±114.9 

23.7 ±102.8 

8.3 ±100.2 

3 

17(a) 

20.2 

9.5 

23.1 ±114.2 

17.1±113.1 

8.3 ±100.9 

1.8±99.9 

3 

17(b) 

S/A 

S/A 

S/A 

S/A 

6.9 ±100.7 

S/A 

S/A 

17(c) 

S/A 

S/A 

S/A 

S/A 

5.3 ±100.3 

S/A 

S/A 

17(d) 

S/A 

S/A 

S/A 

S/A 

-5.0±99.9 

S/A 

S/A 


Table 1: Sky value on the figures in units of counts, counts 5 ^ data multiplied by lO^s. For mock images, the background (= lOOOOe ) was subtracted. LS: low 
symmetricity, see Appendix C. S/A: same as above. N/A: not applicable. 


is accurately known (in a mock image) is the most basic test to 
evaluate the success of a detection algorithm. This comparison is 
done in Section 5.1. Purity is a quantitative measure to assess the 
contamination of false detections that remain in the final result. 
In Section 5.2 purity is defined and used for an analysis of false 
detections in the faintest limits. Finally, in Section 5.3, the faint 
end magnitude dispersion is studied. The statistical significance 
of the imposed requirements on the objects detected as true is 
discussed in Section 5.4. 

While being demonstrative for a “proof-of-concept” paper 
like this one, tests using mock objects or noise can never simu¬ 
late all real circumstances accurately. Thus NoiseChisel is tested 
using only real data without any modeling in a companion paper 
by M. Akhlaghi et al. (2015, in preparation). Its accuracy is fur¬ 
ther evaluated, for example completeness, purity, number counts 
and etc on real data in that work. The details of the real data 
set, necessary modifications to the existing pipeline, tests, and 
an analysis of their results is beyond the scope of this definition 
or “proof-of-concept” paper. 

5.1. Sky Value 

The Sky is the average of undetected regions. Therefore the 
success of a detection algorithm can be assessed with the effect 
of undetected objects on the Sky value. The more successful the 
detection algorithm is, the closer the Sky value will be to the 
background value, which is 10000^“ for all the mock images 
(see Section 2.2). There are two types of “faint” pixels that can 
be left undetected. (1) When the brightest pixels of an object are 
faint, for example. Figure 3(c.2). (2) When the object is bright 
enough to be detected but has wings that penetrate into the noise 
very slowly (for example Figure 3(a.2)). The failure to detect 
both will be imprinted into the average of undetected regions of 
the image. 

Table 1, has all the relevant data for this comparison. Ex¬ 
cept for the last column, all columns in Table 1 show the pixel 
statistics of images shown in the figures of column 1. In these 
statistics, only the pixels within the 200 x 200 pixel box are used 
and not the whole image (see Section 2.1). As a-clipping is still 
extensively used in existing pipelines (see Section A. 3) in Table 
1 the results of the two major a-clipping strategies are also in¬ 
cluded, by convergence and (4,5). In SExtractor, which uses the 
former, the average is used as the Sky value and in the latter, the 


median is used, so we take the same approach here. 

The a-clipping results show that due to the skewness caused 
by galaxies (which do not have a sharp cutoff), none are signifi¬ 
cantly different from the median. Since the final average is used 
in the former, its reported values can be even larger than the ini¬ 
tial median. The correlated noise, which acts like a convolution 
and skews the distribution, causes this difference with the mode 
to be much larger in the real (processed) images of Figures 13 
and 19 (a)-(d) than the mock images. 

In order to check the sensitivity of NoiseChisel, it is applied 
to a set of mock and real images of Figures 12 and 13 with the 
default parameters discussed in Section 3 and Section 4 and fully 
listed in Appendix E. The values for the Sersic index in Eigure 
12 (^ = 0.5,1,4 and 10) were chosen to approximately cover 
the range of observed galaxy profiles (for example, see Lack- 
ner and Gunn 2012). The objects in Eigure 13 are HST/ACS, 
F814W images from the COSMOS survey of what are thought to 
be three z ^ 1.5 star-forming galaxies, so that they are imaged 
in rest-frame UV. The same input images are given to SExtrac¬ 
tor (Eigures 18 and 19 in Appendix B with the configuration file 
shown in Appendix D) and the results are reported here. Some 
studies use the average of undetected pixels when SExtractor’s 
detections have been removed. Therefore in column 6 of Table 1 
the average pixel value when all SExtractor detections have been 
masked to 3,3.is reported. 

The mock images (with a number in the last column) show 
that compared to a-clipping, removing SExtractor’s detected ob¬ 
jects has been more successful in approaching the real back¬ 
ground value (0e~), but the average of the undetected regions 
are still a significant overestimation. The last 4 rows of Table 1, 
which successively mask more of the bright profile, show that 
this is not only due to undetected objects, but also the fainter 
parts of bright objects. The last row shows that when too much 
of the image is masked, the average of non-detected pixels goes 
below the real Sky value. This is reasonable because localized 
high valued noise pixels are systematically removed. 

Over-estimating the Sky value results in a systematic under¬ 
estimation of the total count of the detected objects. In the most 
diffuse and low surface brightness cases, a severe Sky overesti¬ 
mation can even result in a negative total count for the targets. 
Some of the individual pixels of a real object can be below the 
Sky value (Section 3.1.2) but if the Sky value is accurately found, 
the total count of a real detection can never be negative. If it oc- 
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Figure 12: NoiseChisel sensitivity test. In each image there are 25 mock pro¬ 
files, equally spaced (see Figure 3(b.l) and (c.l) for the images prior to adding 
noise for the cases of n = 4 and n = 1). Except for their total magnitude, all 
profiles in each image have the same parameters, = 10 pix, 0 = 45°, q = senq. 
The Sersic index in each input image is different. The profiles on the bottom 
left are the faintest with —9.2 magnitude while those on the top right are the 
brightest with —10.64 magnitude. The second column shows the regions of the 
detected objects masked from the input image. The third column is the inverse 
of the second. 


curs, it can only be due to the overestimation of the Sky. This 
problem can only occur in the existing signal based approach to 
detection where any pixel above the threshold is assumed to be 
the top of a profile. Negative total counts have been reported 
and extensively used in the SDSS survey, for example. However, 
only the problem of not being able to measure the logarithm of 
a negative value (for conversion to a magnitude) was addressed 
(see Lupton et al. 1999). The source of the problem, which is 
an overestimation of the Sky, was not sought. In NoiseChisel, if 
the numerator of any of the S/N measurements (Equations (3)- 
(5)) becomes zero or negative, that detection, clump, or river is 
discarded. 

For the mock galaxies, where the background is accurately 
known, when the area covered by a diffuse source is the greatest 
(Figures 12(a), 13(e) and 17(a) and (b)), NoiseChisel’s detection 
algorithm is 2.3, 2.9, 4.6, and 3.8 times more successful than 
SExtractor. 

For Figures 13(e) and 19(e), Table 1 shows the edge and 
crowded field effects on both SExtractor and NoiseChisel. Com¬ 
pared to the other mock images, NoiseChisel’s Sky value is sys¬ 
tematically higher in this case. This is not due to the edge but 
due to the six large ^ = 4.00 profiles that exist in this image (see 
Figure 1 (b. 1) for the complete pixel coverage of one of those six 



Figure 13: NoiseChisel applied to realistic conditions, (a)-(d) Four real galax¬ 
ies. (e) Six mock profiles with the same parameters of Figure 1(b) and position 
angle of 0 = 45°. Five of the profiles are centered 35 pixels outside the bottom 
and left edges of the image; see Section 2.1. The columns are the same as those 
in Figure 12. The truncation count (where black pixels are defined) in (e.l) is 
half that of Figure 1(b) to emphasize the faint wings of the profiles outside the 
image edge. 

profiles). If — smesh was set to a larger size to find the threshold 
(see Section 4.2.1), the meshes covered by the objects would be 
ignored and NoiseChisel’s final Sky value would be lower. How¬ 
ever, Figure 11(d) shows that even such a small —smesh was not 
enough to completely remove some strong gradients. Therefore 
this higher Sky value can be corrected if there are no gradients 
in the image. Only one value for the full displayed region was 
discussed here. In practice, SExtractor interpolates over such re¬ 
gional values to find a background value for each pixel. Section 
A.6 discusses the process and shows the histogram of the pixel 
by pixel Sky values used. 

5.2. Purity 

Assume that we have made mock profiles in mock noise 
where the number and positions of the mock objects buried in 
the noise are already known. By construction we know there 
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Figure 14: Comparison of purity (false detections) and dispersion in measured magnitude when completeness is 1 for identical and very faint n = 0.5 profiles, (a) and 
(b) 25 profiles are located in the central 200 x 200 region, so any detection outside of it is false, see Section 2.1. Completeness is constrained to 1 for both techniques. 
SExtractor and NoiseChisel have 69 and 7 false detections respectively. The two histograms (c) and (d) show the distribution of magnitude (with zeropoint magnitude 
of 0) for the true and false detections for the respective technique. The dashed line shows the true magnitude of all the profiles. The histograms are derived from four 
runs (100 true detections) on an identical no-noised image but with different simulated noise seeds, (a) and (b) had the most number of false detections for each method. 


is no other source of signal in such an image. Take Nf to be the 
number of false detections, Nt, the the number of true detections, 
Nt, the total number of detections, and Nj the total number of in¬ 
put mock profiles. Therefore Nf^Nt = Nt. The completeness, 
C, and purity, P, of a given detection algorithm with a given set 
of parameters are defined as 





( 6 ) 


Completeness is commonly measured with mock profiles, 
for example, mock point sources or mock galaxies. The results 
from such mock objects are not realistic. This is because, as dis¬ 
cussed throughout this paper, real galaxies display very diverse 
morphologies and do not generically satisfy the clean elliptical 
and radial profiles that can be modeled. Therefore any result 
based on mock objects will overestimate C. When the detection 
algorithm also depends on profiles having simple shapes with a 
clearly defined center and radially decreasing profile, the sys¬ 
tematic bias is exacerbated. Purity, on the other hand, deals with 
false detections which don’t require modeling. 

Without true detections, purity will be 0 since = 0 or 
Nf = Nt. In general, for a given set of parameters in any detec¬ 
tion algorithm, C and P are anti-correlated: as purity increases, 
completeness decreases and vice versa. Therefore neither can 


be studied independently. For example, in SExtractor decreasing 
the threshold in Figure 18 will allow the detection of more of the 
faint mock profiles therefore increasing completeness. Simulta¬ 
neously, however, the number of false detections will increase 
(see Figure 17, for example) and thereby decreasing the purity 
of the output. 

The methodologies of NoiseChisel and SExtractor are fun¬ 
damentally different, therefore comparing their purity should be 
done carefully. To have a reasonably objective measure of purity 
in this case, the only way is to constrain completeness between 
the two. Therefore a mock image that has C = 1 with the default 
parameters of NoiseChisel is created. All the parameters of SEx¬ 
tractor are fixed to Appendix D, except detect_thresh which is 
decreased until C = 1 is achieved in SExtractor too. This strategy 
is chosen based on the advice of SExtractor’s manual to keep the 
sensitivity only related to detect_thresh. The result will not be 
independent of the mock profiles used to define completeness. 
The sharper the profiles, the fewer faint pixels will be needed to 
reach C = 1. Hence an ideal test would be on the fiattest profiles 
or those with n = 0.5, which both SExtractor and NoiseChisel 
show as their lowest completeness (see Figures 12 and 18). 

To get an accurate measure of purity, a large area of the im¬ 
age has to have no data. The mock image is very similar to the 
input image of Figures 12(a.l). Recall that 96% of the area of 
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all the mock images in this paper is blank noise (Section 2.1). 
The only difference is that all 25 profiles have the same total 
— 10.5 magnitude. This is done to decrease the effect of ran¬ 
dom noise positioning. The fixed magnitude of the mock profiles 
is set to be such that with the default parameters of Appendix 
E, NoiseChisel detects all 25 in the image. To make the test 
more robust, the mock noisy image is created four times, such 
that only the noise differs. The magnitude for all 25 x 4 = 100 
profiles was the faintest magnitude that NoiseChisel would give 
C = 1 in all four images. 

S Extractor could not detect all the objects without fiagging 
some of their magnitudes as unreliable. So ignoring the flags, 
the largest threshold that S Extractor also gives C = 1 for all 
thresholds below it and for all four mock images is found to be 
DETECT_THRESH=0 . 6 . The masked results of one of the images 
can be seen in Figure 14 once masked with S Extractor’s detec¬ 
tions and once with NoiseChisel’s. In total, in the four images, 
100 mock galaxies are present and in all four, completeness is 
1. Therefore the result of all four images can be reported in one 
purity value for SExtractor (Ps) and NoiseChisel (Pn). 


25 + 25 + 25 + 25 _ 
79 + 90 + 94 + 87 


Pn = 


25 + 25 + 25 + 25 
27 + 26 + 32 + 27 


= 0.89 


Even though the threshold used in NoiseChisel is so much 
lower than the smallest possible threshold to SExtractor (which 
is the Sky value), this result shows that the purity of NoiseChisel 
is about three times higher in detecting real, faint objects with the 
parameters used (see Appendices D and E). In the last column 
of Table 1 the number of false detections (Nf) in all the mock 
images is reported. Recall that like Figure 14, all mock images 
are also covered by Gaussian noise for 96% their area. With 
exactly the same input parameters, the reported Nf values agree 
well with those reported in this test. 

As discussed in Sections 3.1.5 and 3.2.1, NoiseChisel is very 
robust in dealing with correlated noise that processed real images 
also have. However, SExtractor and all existing signal-based de¬ 
tection methods in general treat the pixels above the threshold 
as the top of a known profile; therefore they are not immune to 
the effects of correlated noise, resulting in lower purity when 
applied to a real processed image. 


5.3. Magnitude Dispersion 

All the profiles in Figure 14 are identical, however. Figure 
14(a) displays a large dispersion in the shapes and areas in the 
detections. This results in the roughly 2 magnitude rage that is 
visible in SExtractor’s true detections (see Figure 14(c)). It is 
very important to note that the bimodality that is observed in 
Figure 14(c) (ignoring color) is due to the fact that all mock pro¬ 
files were identical. In a real image, the faintest true detections 
can have a variety of inherent profiles. Therefore the distribu¬ 
tion will be unimodal with negative skew. The net result is a 
common plot in existing catalogs (for example see Figure 10 in 
Illingworth et al. (2013)). On the contrary, NoiseChisel’s true de¬ 
tections all have roughly the same area in Figure 14(b) and thus 
the range in the measured magnitude for true detections in Fig¬ 
ure 14(d) is about half that of SExtractor. The two NoiseChisel 


detections with magnitudes brighter than — 11 are the result of 
four of the profiles being detected as two detections. 


5.4. Imposed Restrictions on the Signal 

The proposed algorithm imposes some restrictions on the ob¬ 
jects it will detect as true. Therefore it is not “absolutely” noise 
based or independent of any signal parameter. Here we will dis¬ 
cuss its methodological limits. With the chosen parameters in 
this paper (listed in Appendix E) the target has to satisfy the fol¬ 
lowing conditions. 

1. To be detected, an object has to have an area equal to or 
larger than Figure 6(c), above the initial quantile threshold 
(which is below the Sky value; see Section 3.1.2) on the 
convolved image. 

2. To be classified as a true detection, it has to contain at least 
one connected component larger than 15 pixels above the 
Sky threshold with a sufficiently large S/N specified from 
the ambient noise (see Section 3.1.5) on the actual image. 

As the methodological limits or requirements imposed on the 
target by a detection technique decrease, the detection algorithm 
becomes more successful because it becomes more generic. The 
“success” of a detection algorithm can be defined as how accu¬ 
rate its Sky measurement can approximate a known background, 
the purity of its result and the scatter in its photometry as dis¬ 
cussed above. It was shown that because of far fewer constraints 
on the target objects, this technique was significantly more “suc¬ 
cessful” compared to SExtractor (with the parameters of Ap¬ 
pendix D). 

In astronomical data analysis, an S/N above five is usually 
considered as a “solid” result, any measurement with a lower S/N 
is usually included in an analysis with caution (McLean 2008, 
Section 1.5.2). Using ideal simulated noise we showed that ob¬ 
jects and clumps with an S/N as low as 3.61 (Figure 7) can be 
accurately (0.99%) detected with this technique using this set of 
parameters. Note that the Sky value used in those plots was a 
first (under-estimated) approximation, therefore the counts and 
hence the S/Ns are overestimated. Hence the true limiting S/N is 
slightly less than what is reported in those plots. 

Through setting fewer constraints (if a similar purity level 
can be achieved) it would be possible to decrease the S/N of true 
detections even further. The extreme case would be to detect ob¬ 
jects with an S/N of unity. However, the statistical significance 
of any resulting measurement (for example, the effect of that de¬ 
tection on the measured Sky, its total count, its central position, 
or any other measurements on the object) will correspondingly 
decrease to values that are no longer significant in any scien¬ 
tific/statistical measurement and analysis. 

6. DISCUSSION 

A new method for detecting extremely faint objects and fainter 
parts of brighter objects in noise is introduced. Unlike the exist¬ 
ing approach to detection where the signal is the basis, this ap¬ 
proach is based on operations and calculations on the noise with 
insignificant requirements for the signal. This method is also 
non-parametric, in other words, it does not involve a functional 
regression analysis. It is therefore ideal for the study of nebu¬ 
lous/amorphous objects that are heavily immersed in noise, for 
example, galaxies. 
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Galaxies can have rich dynamic histories involving internal 
processes as well as external interactions with the halo or other 
galaxies. This creates very complicated, diverse morphologies. 
The fact that we cannot image them in their full 3D glory further 
complicates their final observed image. Therefore, imposing any 
a priori model on their shape or radial light profile during the de¬ 
tection process is a self-inflicted systematic bias in their study. 
Modeling the light profiles of galaxies is vital to our understand¬ 
ing of the galaxies and testing our models, but it should be done 
after detection is complete. Mixing detection and modeling will 
only bias both. 

One major argument in favor of elliptical-based detection 
and photometry techniques like the Kron and Petrosian radii 
as opposed to isophotal methods (including NoiseChisel) is that 
they are impervious to (1 dimming and K-correction and 
for any profile, they are “well defined.” By well defined it is 
meant that if the same galaxy (with the same morphology) exists 
in multiple redshifts, then such methods will yield the same frac¬ 
tion of total light for all the redshifts. Ideally, when the apertures 
used extend to infinity (see Graham and Driver 2005), and the 
galaxies have no morphological evolution, the result would be 
independent of the noise and thus the threshold used. However, 
this is not necessarily the case in the real world. 

The interpretation of signal-based detection results depends 
on the morphology of the galaxies under study. To demonstrate 
this point, consider the simple mean in a ID distribution as an 
analogy. Note that the Kron radius is a weighted mean. The 
mean is a well defined statistic for most distributions because it 
gives a unique and unambiguous value. Regardless of whether 
the underlying distribution is a Gaussian or a log-normal dis¬ 
tribution, for example, an unambiguous mean can be defined, 
but for interpreting the mean, the distributions should not be ig¬ 
nored; otherwise the analysis will be biased when such different 
distributions are involved. 

The same applies to the Kron or Petrosian radii. For any 
given surface brightness (morphology), a unique and unambigu¬ 
ous radius can be found, rendering them well defined. How¬ 
ever, when it comes to interpreting what that radius physically 
means (in kiloparsecs or fraction of total light, for example), the 
interpretation will be biased when the morphology (spatial and 
count distribution of pixels) differs between the sample of galax¬ 
ies. Therefore, comparing Kron magnitudes for the galaxies of 
a general population can be problematic. For example, all the 
galaxies in a (noisy) survey or image which can simultaneously 
contain irregulars, spirals, and elliptical galaxies. 

The Kron or Petrosian radii satisfy the redshift independence 
in studies of galaxies that are already known to have approx¬ 
imately the same morphology; in other words, a morphologi¬ 
cally selected sample where the sample has been detected in¬ 
dependently of the shape of the galaxies. For example, unlike 
the references in Section 1, galaxies might indeed have simi¬ 
lar morphologies across redshifts (for example, see Lee et al. 
2013). However, this claim, or any other claims, including those 
in Section 1, about the morphology of detected sources can only 
be tested when the detection technique does not depend on the 
morphology as an assumption. In other words, if the detection 
technique depends on galaxies having a similar morphology at 
similar or various redshifts, it should not come as a surprise if 
galaxies with similar morphologies are detected. 

Isophotal detection methods like NoiseChisel do have the 
limitation that they use a different threshold for objects at differ¬ 


ent redshifts, or {Iz)^ dimming. However, it should be con¬ 
sidered that the Sky value and its standard deviation in the input 
image are also completely independent of such issues. These 
vital image statistics, which propagate in all subsequent mea¬ 
surements, only depend on the arriving photons in each pixel 
regardless of the physical origin of the source. Therefore im¬ 
posing such high-level (related to physical interpretation and not 
the data) constraints on low-level measurements (to do with the 
data, regardless of later physical interpretations) is a source of 
systematic bias. This noise-based approach to detection was de¬ 
signed to create an accurate low-level measurement tool with as 
few assumptions as possible so high-level interpretations can be 
less biased. 

High-level interpretations of the outputs of low-level mea¬ 
surements can easily be done after the low-level procedure is 
complete. For example, to account for (1 dimming, an in¬ 
dividual threshold can be defined based on the measured redshift 
(a high-level product) for each object in a study. In the 1970s, 
when the Kron and Petrosian radii were defined, the processing 
power to achieve this level of customization was not available. 
Therefore such high-level constraints had to be imposed on the 
techniques that originated decades ago. However with the very 
fast computers cheaply available today, there is no more need for 
such self-inflicted biases. 

Through the particular definition of threshold (in Figure 2, 
Section 3.1.2 and 4.2.1) the cycle of the Sky estimate depending 
on detection and detection depending on Sky are significantly 
weakened. In the proposed technique, the detection threshold 
and initial detection process are defined independently of the 
Sky value. An independent initial approximation of the Sky 
is found by averaging the initial undetected regions and is only 
used in classifying and removing false initial detections (see the 
flowchart of Figure 2). 

Due to all the novel, low-level and non-parametric methods 
employed in this noise based detection, fainter profiles can be de¬ 
tected more successfully than the existing signal-based detection 
technique (at least with the parameters used here). The purity 
measurements of Section 5.2 show how successful this algorithm 
and its implementation have been in removing false detections, 
while detecting the very faint true profiles and the fainter parts of 
large profiles. Section 5.3 then showed the significantly smaller 
scatter in the photometry of those very diffuse and faint galaxies. 

A new approach to segment a detected region into objects 
that are heavily immersed in noise is also developed. In this 
approach, the true clumps are first found using the global prop¬ 
erties of the noise, not based directly on a user input value or 
the particular object (compare Section 3.2.1 with Section B.2). 
By expanding the clumps over the light profile of the detected 
region, possible objects are found and the detected region is seg¬ 
mented. 

One of the design principles of NoiseChisel was to decrease 
the dependency of the output on user input. To ensure that the re¬ 
sults are more dependent on the image than the user’s subjective 
experience. To achieve this, the true/false criteria were defined 
based on the image and not an absolute hand input value (see 
Section 3.1.5 and Section 3.2.1). This approach accounts for 
local variations in the noise properties or correlated noise that 
becomes a significant issue in processed data products. Note 
that by using relative measures the dependence on user input has 
significantly decreased, but not disappeared. 

This noise based approach to detection allows for the input 
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parameters to require minimal customization between images of 
different instruments and in different processing stages. For ex¬ 
ample, in this paper, all the outputs for mock images with ideal 
noise, raw ground-based images with real noise and HST/ACS 
drizzled space-based images with correlated noise were processed 
with nearly identical input parameters. 

NoiseChisel is designed with efficiency in mind. Therefore 
it is also very fast in processing a large image. As an exam¬ 
ple the image of Figure 11 was processed in 4.43 seconds using 
a 3.06GHz, 4 physical core CPU on a desktop computer. Until 
now the design and applicability of the concepts introduced here 
was the primary concern, but in time more efficient algorithms 
and methods might be found to increase its efficiency. For ex¬ 
ample including GPU functionality can significantly decrease the 
running time. The detection and segmentation techniques intro¬ 
duced through NoiseChisel are easily expandable to the detection 
of any signal in noise, for example, to 3D data in radio studies 
similar to clumpfind (Williams et al. 1994) and ID data like spec¬ 
troscopic data. 

The technique proposed here for detection and segmentation 
along with the ancillary NoiseChisel program provides a funda¬ 
mentally fresh and new approach to astronomical data analysis. 
The limited preliminary tests on mock profiles and noise in this 
paper show that it significantly increases the detection ability, 
with a reasonable purity and very small photometric scatter, to 
detect the fainter parts of bright galaxies with any shape or to 
find small faint objects that will be undetected with the existing 
signal based approach to detection. Both of these new abilities 
can play a major role in most branches of galaxy evolution or a 
astronomy research as whole (some applications are reviewed in 
Section 1). Along with the new instruments such as James Webb 
Space Telescope, the Large Synoptic Survey Telescope (LSST), 
and the Thirty Meter Telescope, which will be available in the 
next few years, this new approach to detection will play a very 
important role in current and future astronomical discoveries. 
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A. EXISTING SKY CALCULATION METHODS 

In existing algorithms for detection, the first step of detec¬ 
tion, namely setting a threshold, is based on the Sky value (see 
Appendix B). To approximate the Sky value prior to detection, 
researchers use multiple techniques, such as the mode, a-clipping, 
etc. A critical analysis of the most common techniques for ap¬ 
proximating the Sky value is given here. 

A.i. Local Mode 

In a symmetric random distribution, for example, pure noise 
in an astronomical image, the mean, median, and mode are 
equal within their statistical errors. With the addition of signal 
(astronomical objects, yellow histogram in Eigure 1), a skewed 
distribution is formed (red histogram in Eigure 1). The mean is 
the most highly dependent on the skewness since it shifts toward 
the positive much faster. The median is less affected since it will 
be smaller than the mean for a given positively skewed distribu¬ 
tion. Einally, the mode of the distribution is the least affected 
statistic of the three. This argument has thus led many to take 
the mode of a noisy image as the Sky value. 

There are several approaches to finding the mode of the im¬ 
age: Bijaoui (1980) finds a functional form to the image his¬ 
togram in the vicinity of the mode using Bayes’s theorem. Kron 
(1980) finds the mode of the pixel distribution in a 50pixel x 
50pixel (20" x 20") box about the peak of each object based on 
the mean of the 7 bins in the histogram with the largest num¬ 
ber of pixels and considered that as the Sky value. Beard et al. 
(1990) use a more elaborate approach in the search to find the 
mode. The histogram is first smoothed with a moving box filter. 
The highest peaks are found and used in a cubic fit to find the 
mode which they consider to be the Sky. 

DAOPHOT (Stetson 1987) also uses a very similar method: 
a circular annulus with an inner radius several times the stel¬ 
lar FWHM is taken from a nearby part of the image, sufficiently 
close to the target and the mode of the pixels in that region is 
considered as the Sky value. The mode is approximated based 
on the following relation between the mean and median: 3 x 
median — 2 x mean. This particular relation between the mode, 
median and mean, can only exist when one assumes a particular 
pixel count probability distribution function (PDF) or histogram 
for the image pixels. 

SExtractor can be considered a hybrid, first relying on a- 
clipping, and then on finding the mode to find the Sky value. 
The former will be discussed in Section A. 3. In the latter, like 
DAOPHOT, SExtractor uses the following relation to find the 
mode: 2.5 x median — 1.5 x mean. In a recent re-analysis of the 
SDSS stripe 82 (the SDSS deep field), Jiang et al. (2014) use a 
similar but more accurate hybrid method. The main difference is 

The photon noise actually comes from a Poisson distribution which is not 
symmetric especially for very low mean (k) values, but with the very high 
mean values (k — B— 10000e“ from Figure 1), the Poisson distribution can 
safely be considered approximately equal to a Gaussian distribution with /j = 
k and a = Vk . Images taken from the HST commonly have much smaller 
Sky values. In such cases, the main source of noise is noise produced by 
the instrumentation, for example, read-out noise which can out-weight the 
non-symmetric low-^ Poisson noise. 
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that they remove SExtractor detections first—with detect_thresh=2 
and DETECT_MINAREA=4 (see Section B.l). If the mean is not 
smaller than the median, the mode is assumed to be found based 
on a similar relation to DAO PHOT between the mean and mode. 

Therefore existing methods to use the mode as a proxy for 
the Sky either employ the parametric approximations of Bijaoui 
(1980), DAOPHOT, SExtractor, and Jiang et al. (2014) based on 
an assumed function for the pixel count PDF or a non-parametric 
one (Kron 1980; Beard et al. 1990). It is clear from Figure 1 
that the histogram, which can be considered a binned PDF can 
take any form of skewness depending on the type and position¬ 
ing of astronomical data in the image. Therefore by assuming 
a fixed generic functional form for all possible PDFs, significant 
systematic errors will be induced. 

The existing non-parametric methods reviewed here rely on 
the image histogram. The results from a histogram depend on 
the bin-width and the bin positioning, especially in the vicinity 
of the mode. For example note how the peaks near the mode vary 
in Figure 20 only due to bin positioning. The method of Kron 
(1980) for finding the mode (explained above) for the histogram 
in Figure l(b.l), yields 4.69^“ ± 16.02^“, where the histogram 
is made from 200x200 pixels; recall that Kron (1980) used a 
50 x50 pixels box, and the area here with no signal is much 
larger than the object. Different bin widths would change this 
result, but since no generic distribution can be assumed for all 
images, any generic bin width choice would ultimately be arbi¬ 
trary. A very accurate, non-parametric method for finding the 
mode of a distribution that does not rely on the histogram, but 
on the cumulative frequency plot of the pixels is introduced in 
Appendix C. 

Regardless of how the mode is found, such attempts at using 
the mode of a distribution as a proxy for the Sky fail to consider 
the fact that the mode of the distribution also shifts significantly 
from the actual image background depending on the data that is 
embedded in the noise. Figure 1 shows that like the mean and 
median, the mode is a biased estimator for the Sky value due to 
the data. As the fainter parts of large objects or a large number 
of faint objects cover a larger fraction of the image, the mode 
of the distribution shifts to the positive. Figure 1 shows how 
this shift increases, with mode values of 9.90^“ (symmetricity 
of 1.26, see Appendix C), 9.50^“ (0.48) and, 156.10^“ (0.44), 
respectively. Recall that the background value is 0e~ in all three. 
See Section 5.1 for further discussion on the mode as compared 
to the true Sky. 

Reversing the argument above (that data cause a difference 
between the mode and median), we conclude that if the mode 
and median are approximately equal, there is no significant con¬ 
tribution of signal or data. In fact this argument is very important 
for NoiseChisel (see Section 4.2.1). Note that this argument is 
only possible when the mode is found independently from the 
median as in Appendix C. 

Based on the idea in the previous paragraph we have also 
created SubtractSky, which is also distributed as part of the GNU 
Astronomy Utilities. It can be used for cases where only Sky 
subtraction is desired. Note that the argument above “detects” 
signal through its effect on the distance of the mode and median, 
but only based on pixel values and independently of the signal’s 
spatial distribution. See the manual for more information; its 
operation is very similar to that explained in Section 4 with small 
modifications. 


A.2. Removing Detections 

If the detection algorithm is independent of noise, this method 
of finding the Sky is the most accurate (see Section 2.2). How¬ 
ever, the existing detection algorithms depend on knowing the 
Sky before they run (see Section B.l). Even so, this method is 
used by some researchers to find the Sky value, for example, 
Stoughton et al. (2002) used this technique over a grid to es¬ 
timate the Sky in the SDSS survey. Bright objects, which are 
defined as those with at least one pixel above 200a, are first sub¬ 
tracted from the image. They are all assumed to have power-law 
wings which are modeled and subtracted from the image and 
finally a median filter is applied to the image in a box of side 
100" « 252 pixels. Increasing the box size will be a serious 
burden on the computational process of the pipeline and there 
will always be objects that are, collectively or individually, large 
enough to bias the result. 

Mandelbaum et al. (2005) reported a systematic decrease in 
the number of faint galaxies near < 90" of the center of bright 
objects. It was also observed that this method underestimates 
the total count of the brighter galaxies (for example, Lauer et al. 
2007). The issue was finally addressed in Aihara et al. (2011). 
The major solution was to change the threshold to detect bright 
sources to 51a. Detected objects were deblended, and a linear 
combination of the best exponential and de Vaucouleurs models 
were subtracted. 

While the decreased threshold and new fitting functions do 
increase the number of objects that are masked, they still fail in 
the following respects. (1) Subtraction is done based on para¬ 
metric model fitting which can potentially be severely wrong in 
the fainter outer parts of the galaxies that are the primary source 
of bias. Such cases are when the galaxy cannot be character¬ 
ized by a simple ellipse or when it is on the edge of the image 
(see Section B.4). Other cases are when the galaxy might not 
be a simple ^ = 1 or ^ = 4 profile. (2) The bright and extended 
objects that do not have a pixel above 51a are still a cause of 
systematic bias. For example, none of the bright profiles in Fig¬ 
ure l((c), untruncated) have a pixel above that extremely high 
threshold. The brightest pixel in Figure l(c.l) is only 39.67a 
above the background. 

A.3. a-clipping 

When the outliers of a data set are sufficiently distinct from 
the majority of the (noisy) distribution, a-clipping can be very 
useful in removing their effect. Assuming m is the distribution 
median, any data point lying beyond m^ac can be clipped (re¬ 
moved) and this process can be repeated on the smaller data set 
n times. The parameters to define a-clipping can thus be written 
as Such a-clipping is only useful if the objects (outliers) 

have a very high S/N with very sharp boundaries, for example 
cosmic rays (see Figure 15). 

Some prominent uses of this approach in finding the Sky 
value in an image are: the HST image processing pipeline (Gon- 
zaga et al. 2012) with (4,5) and SDSS Stripe 82 (Annis et al. 
2014) with (3,5) and again (Jiang et al. 2014) with (_^'^,20). As 
its primary tool, SExtractor also relies on this approach though 
the termination criterion is not a fixed number of times, it is the 
convergence of a. Asymmetric clipping was proposed by Rat- 
natunga and Newell (1984). Gonzaga et al. (2012) obtains the 

Jiang et al. (2014) do not to mention the first parameter. 
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Figure 15: Effect of cosmic rays on the histogram in a small postage stamp 
(100 X 100 pixels) of a raw Sky subtracted HST/ACS (fit. fits) image. The 
histogram is vertically truncated and the color range in the displayed image has 
the same range as the histogram. 84 cosmic ray pixels have values above I50e~ 
and are not shown in the histogram. The outlying cosmic rays can easily be 
distinguished in the histogram due to their sharp boundaries in contrast to the 
fading boundaries of galaxies, for example. Figure 1(b). 

value globally and uses it for the whole image; others use this 
technique to find the local Sky value and use interpolation over 
the whole image. The five vertical lines in all three examples of 
Figure 1 show the position of each iteration’s a-clipped median 
with (4,5) on the mock images. Due to their extreme proximity, 
they might not be resolved. 

Once a converges, SExtractor uses the mean as the Sky value. 
For the displayed area of the three examples of Figure 1, the 
mean, when SExtractor’s convergence is achieved, is, respec¬ 
tively, 11.1±103.3^-, 23.1±114.2^-, and 236.9 ± 215.6^“. 
Prior to a-clipping, the median of each image was: 11.1^“, 
20.2^“, and 199.3^“. Hence, a simple median (no a-clipping) 
was a better approximation to the true background value (0 ± 
100^“). The medians of the (4,5) a-clipping technique are shown 
on Figure 1 . It is clear that because the profiles of astronomical 
objects go slowly into the noise, the correction provided by a- 
clipping is very insignificant for the mean and median. For the 
standard deviation, however, a-clipping is very useful. The ini¬ 
tial standard deviation of these three examples was 167.30^“, 
455.90^“ and 341.05^“ (see Section 5.1 for a more complete 
comparison). 

Figure 1 shows how through their very gradual penetration 
into the noise, stars (the PSF), galaxies, and undetected objects in 
astronomical images are very different from cosmic rays. There¬ 
fore a-clipping will not suffice when such objects are present. In 
other words, no matter how much the brightest pixels of objects 
in an image are clipped, their faint wings penetrate very deep 
into the noise. Therefore when astronomical objects are present 
in the image, a generic a-clip is a very crude overestimation of 
the Sky value. 

A.4. Radius of Flat Profile 

Another possible solution to finding the Sky was presented 
recently as a proxy for the local Sky background behind an ob¬ 
ject. Barden et al. (2012) has proposed finding the a-clipped 


mean value on elliptical annuli after removing all known source 
contributions, using SExtractor. As long as the a-clipped mean 
is very high above the noise, the slope of the mean pixel value 
on an annulus as a function of radius will be negative. However, 
as it approaches the noise, due to the scatter caused by noise, 
the slope will become repeatedly positive and negative. They 
consider the position of the second positive to be the annulus on 
which they find Sky value for each object. 

The drawbacks to this method are that it assumes that the 
profiles can be represented as an ellipse with a well-defined cen¬ 
ter. It is also based on the Kron radius (Section B.4) which can 
underestimate the full extent of an object when detect_thresh>1 
(see Section B.1.2 for a complete discussion). The most impor¬ 
tant issue with this approach is the area (or number of pixels) 
which is used to estimate the Sky value. As Stetson (1987) noted, 
to reach a statistically acceptable precision, the area used to find 
the Sky value has to be sufficiently larger than the object itself. 
However, the average value on an annulus is based on a smaller 
number of pixels than that spanning the parent object especially 
when all neighboring objects are removed (see their Figure 5). 
Therefore the random scatter due to the small number of pixels 
used will limit the ability to measure the true Sky level with ac¬ 
curacy. This is further exacerbated by the fact that the growth 
of the elliptical annuli is halted based on scatter in the measured 
Sky of each. 

A.5. Median Value of Dithered Images 

Median values of dithered images can only be used in images 
prior to co-adding. This method exploits the dithering, where the 
field of view of the telescope is shifted slightly between differ¬ 
ent exposures (for example in Kajisawa et al. 2011). In order 
to find the Sky value between several temporally close observa¬ 
tions, the median pixel value of those observations is taken as 
the Sky value on each pixel. 

The weaknesses of this method are: (1) the Sky value might 
change between separate observations, particularly in infrared 
imaging. (2) The main objects in the image have to be far smaller 
than the dithering length. (3) The field should not be crowded. 
(4) It assumes perfect bias, dark, and fiat fielding since it depends 
on operations that are done prior to them. 

A. 6. Background Interpolation 

SExtractor and most other algorithms do not find one Sky 
value for the whole image. Instead, the Sky values are found 
on a mesh grid. In SExtractor the sizes of the meshes are set 
to BACK_siZE pixels. It then uses bicubic-spline interpolation to 
find a the Sky value for each pixel. Figures 16(a)-(c) show SEx¬ 
tractor’s Sky value calculated for each pixel for the three mock 
images of Figure 1 along with the images in Figures 18 and 19. 
Recall that the background value for all mock images is zero. It 
is clear that as the area covered by the faint parts of mock galax¬ 
ies has increased from (a) to (c) the background values calcu¬ 
lated by SExtractor’s a-clipping and mode approximation have 
been shifted very strongly to the positive. In this figure the prob¬ 
lems with using model based interpolation techniques are also 
evident. Note how the corners of the images have extremely 
high and low values compared to the centers of the images in 
Figure 16 (a)-(c) and (1). Note that the other cases of Figure 16 
do not include image corners. 
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Figure 16: Sky value check images and their histograms, produced by SExtractor with the input parameters of Appendix D applied to the images in some of the figures 
in this paper. First row (a)-(c): the three mock images of Figure 1. Only for these three, the displayed 200 x 200 region was input to SExtractor. Second row (d)-(g): 
the four mock images of Figure 18. Third row: (h)-(l) the four real and one mock image of Figure 19. Failure to interpolate accurately on the corners of the images are 
clearly visible in (a)-(c) and (1) where the image edge can be seen. Real images are in units of counts s“^ so that the horizontal axis of their histograms are multiplied 
by 10^5. The true background (B) value of all mock images is subtracted in the histograms. 


In SExtractor, through the parameter back_size (see Appendix 
D), a user can set a smaller grid on the image. SExtractor then 
finds the Sky in those sub grids, removes the outliers with me¬ 
dian filtering, and interpolates over them. Therefore for each par¬ 
ticular image of a particular target, a good choice of back_size 
will slightly correct the very large overestimations discussed here. 
However, for any back_size, there might be objects in the image 
that cause a similar or even larger bias. Therefore, unless the best 
BACK_siZE is found for each particular input image—depending 
on the objects and their position in the input image—there will 
be a bias. Such customized application is not practical because 
of the sheer number of galaxies studied in most papers. If the 
user chooses a custom background for each image of a large 
set of targets, comparison between the different images cannot 
be accurate because of the different statistical properties of the 
background. Therefore, such significant systematic overestima¬ 
tions will plague any astrophysical analysis or hypothesis testing 
that is based on them. 

B. EXISTING DETECTION AND 
SEGMENTATION/DEBLENDING METHODS 

The detection technique used in astronomy up until now can 
generally be classified as a signal-based approach with various 
implementations. It is impractical to review all implementations 
and their minor differences here. Thus in this appendix we use 
SExtractor 2.19.5 (Bertin and Arnouts 1996) while mentioning 
the important alternative implementations during the steps. SEx¬ 
tractor was chosen here because it can easily be installed and has 
a relatively complete manual. Because of this it is by far the 


most commonly used tool by most authors in generating large 
catalogs or source extraction of known targets and is thus widely 
recognized by the community. Its techniques are also conceptu¬ 
ally very similar and inherited from other packages, for example, 
DAOPHOT (Stetson 1987), clumpfind (Williams et al. 1994) and 
Kron(1980). 

To extract sources from a noisy image, SExtractor first needs 
to find the noise characteristics or the Sky value and its error for 
every pixel. A comprehensive review of the common methods 
for finding the Sky value are thus provided in Appendix A. Based 
on the noise, a threshold is applied and regions above noise are 
chosen as detections (Section B.l). The detections are then de¬ 
fended (Section B.2). The deblended detections are extended 
beyond the threshold using the Kron radius concept that is crit¬ 
ically analyzed in Section B.4. In Section B.3 a discrete hierar¬ 
chical Markov image modeling method for simultaneous detec¬ 
tion and segmentation of objects in an image is discussed. Since 
it also requires growing the objects afterwards, it is discussed 
before Section B.4. The images are outputs of SExtractor. The 
interpretations of how it has operated are based on its manual or 
corresponding paper. The configuration file used to run SExtrac¬ 
tor is shown in Appendix D. 

B.l. Detection 

After finding the Sky value, SExtractor uses it to define a 
threshold (see Section B. 1 .2). The threshold is applied to a smoothed 
image (see Section B.1.1). Of the various regions found above 
the threshold, those with an area smaller than a user defined 
value are excluded as a false detection (Section B.l. 3) and the 
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remaining objects (true detections) are deblended (Section B.2). 
Using the Kron (1980) radius concept, the deblended detections 
are grown on a fixed ellipse to include pixels below the threshold 
(see Section B.4). 

BAA. Convolution 

A convolution kernel comparable with the PSF has been the 
best choice for this purpose. A qualitative argument in support 
of this choice is that it is the widest convolution kernel that can 
be used such that no object is smoothed out since no object can 
be smaller than the image PSF (see Bijaoui and Dantel 1970, for 
a detailed quantitative discussion). Based on this argument, for 
the example runs of SExtractor on the mock images, the same 
PSF (Moffat function with P = 4.76 and FWHM=3 pixels) that 
was used to create the input image was used. In SExtractor the 
convolution kernel is specified by the parameter filter_name. 

B.J.2. Threshold 

In order to detect objects above the noise, a threshold needs 
to be applied to the convolved image. Only pixels with val¬ 
ues above this threshold are considered for object detection. A 
higher threshold will decrease the probability of false detections, 
because fewer noise pixels will have a count above it. However, 
it will miss fainter objects and the fainter parts of brighter ob¬ 
jects. 

In S Extractor, the threshold is defined in terms of the approx¬ 
imated Sky (s, see Appendix A) and its standard deviation (a): 
5'+DETECT_THRESHxa which are measured on the original input 
image (prior to convolution). The values set for detect_thresh 
are different in various studies. The recently released 3D-HST 
WFC3-selected Photometric Catalogs in the five CANDELS/3D- 
HST Fields (Skelton et al. 2014) set this parameter to 1.8. In the 
K^-selected Catalog in the COSMOS/ULTRAVISTA field (Muzzin 
et al. 2013), it is set to 1.7. As a final example, in the Sub- 
aru/MOIRCS deep survey of the GOODS-N field (Kajisawa et al. 
2011), it is set to 1.3. 

In order to go deeper in the noise, Rix et al. (2004) proposed 
a two-stage detection process. In their technique, objects are de¬ 
tected in two runs (values from Rix et al. 2004): (1) a “cold” 
run where detect_thresh is set to a very high value of 2.3, and 
only the brightest regions of the brightest objects are detected 
and (2) a “hot” run where it is set to a lower value of 1.65 to 
detect the fainter objects. The two resulting catalogs are then 
compared and “hot” objects that lie over the “cold” ones are re¬ 
moved from the final catalog. In their “hot” run, such techniques 
use the lowest possible threshold they can consider. Working on 
deeper images, Leauthaud et al. (2007) used 2.2 and 1 for their 
cold/hot detection thresholds, respectively. 

The threshold is calculated based on the original image, but 
it is applied to the convolved image. Therefore it no longer has 
the statistical properties that any multiple of the standard devi¬ 
ation is expected to have, for example, that ^ 95% of the noise 
pixels have a value below 2a. This occurs because of the de¬ 
creased dynamic range after convolution (see Section 3.1.1 and 
Figure 4). In Figure 4, la, which is the standard deviation of 
the un-convolved image, corresponds to 100^“. However, when 
applied to the convolved image of Figure 4(b) only 44 (of the 
40,000) pixels are above this threshold. Note that the kernel 
used in Figure 4(b) is the same PSF that was used to make the 


mock image. In Figure 4(c)-(e), 100^“ far exceeds the brightest 
convolved pixel. 

Based on this idea, Galametz et al. (2013) adopt a different 
approach to the cold and hot detection technique. detect_thresh 
is approximately similar between the two (0.75 and 0.7), but the 
convolution kernels (see Section B.1.1) used are significantly 
different such that the cold run has a much wider kernel (see 
their Appendix A for the respective SExtractor parameters). 

A different approach to finding the threshold of an image, 
called the “False-Discovery Rate,” was introduced in astronomy 
by Miller et al. (2001). In this approach the average fraction of 
false pixel detections to true pixel detections is held fixed when 
assuming a fixed null hypothesis (background and its noise). 
It was implemented as a detection technique in Hopkins et al. 
(2002). All the pixels in the region must be sorted and a p value 
has to be calculated from an approximated mean and standard 
deviation. Thus this thresholding technique also needs the Sky 
and its error prior to the actual detection process. 

B.L3. Detection: True/False 

Once the threshold of Section B. 1.2 is applied to the smoothed 
image of Section B.1.1, the image is divided into two groups of 
pixels: those that are above and below the threshold. The non¬ 
white pixels in Figure 17 (a.l)-(d.l) and column 2 in Figures 
18 and 19 show this image with a threshold applied (they are all 
SExtractor’s segmentation maps). In Figure 17, SExtractor is run 
on the mock image of Figure l(b.2) with the configuration file of 
Appendix D. All of SExtractor’s parameters in all four cases are 
fixed to those in Appendix D except for detect_thresh, which is 
set to 2, 1, 0.5, and 0.1, respectively. The first observation after 
examining Figure 17 is that as the detection threshold decreases, 
the number of false detections increases: 0, 0, 319 and 11296, 
(in the 1000 x 1000 images, see Section 2.1) respectively. 

The major tool available in SExtractor to define and thus re¬ 
move such false detections is detect_minarea. With this param¬ 
eter the user can specify the minimum area that a true detection 
should have. As the threshold decreases, more and more pixels 
become connected. As an example when detect_thresh=0.1 is 
fixed but DETECT_MINAREA is Set to 10, 25, 50, 75, and 100 pixels, 
8989, 3880, 966, 350 and 260 false detections are found. The 
larger detect_minarea will also result in not being able to de¬ 
tect many compact bright objects that will also be present in a 
real image. The SExtractor manual therefore correctly suggests 
setting DETECT_MINAREA to small values (1-5) so that the convo¬ 
lution kernel and detect_thresh define the sensitivity. 

Figure 18 shows the outputs of SExtractor, configured with 
Appendix D, for extremely faint mock galaxies; see Sections 3.1 
and 5.1. SExtractor’s results on a bright and large profile are 
already analyzed in Figure 17. With detect_thresh=1, only the 
top few pixels of the brightest and sharpest profiles have been 
successfully detected. Note that as discussed in Section B.1.2, 
most surveys use much larger thresholds. 

When a large S/N for detections is required for accurate fit¬ 
tings, for example, in photometric redshifts or SED fittings, some 
researchers might opt for discarding such faint detections. How¬ 
ever, their detection in the image is nevertheless very important 
because if they are not detected and removed, they will cause 
an overestimation in the Sky value, hence systematically under¬ 
estimating the total count of the brighter, high S/N, objects (see 
Section 5.1). 
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Figure 17: SExtractor results of various thresholds with configuration parameters from Appendix D except for DETECT_THRESH. The input image only has the one 
profile of Figure 1(b) (see Section 2.1). The threshold values for DETECT_THRESH are (a) 2, (b) 1, (c) 0.5, (d) 0.1. For each threshold there are two images. Left: the 
segmentation map. Right: the elliptical region within three times the Kron radius (r^^, see Section B.4) is masked (white) from the original image. Notice the darker 
profile regions surrounding the elliptical mask that have not been included in the 3,3.ryt elliptical aperture. The smaller elliptical border marks the aperture containing 
~ 90% of the counts in the profile. The larger elliptical border shows the region containing ~ 95% of the counts. The segmentation maps in (c.l) and (d.l) have 
different color codes for the main object (gray) and any other detections (black). The marked regions show how connectivity is not a criteria in SExtractor’s deblending 
algorithm. 


B.2. Segmentation/Deblending 

Segmentation and deblending are defined in Section 3.2. An 
interesting deblending algorithm is proposed by Lupton (2015) 
and is employed in the SDSS, Subaru Telescope Hyper Suprime- 
Cam, and LSST pipelines. In one image (one color) real peaks 
are defined as those that are at least SCsky above the saddle point 
or valley that separates it from neighboring peaks in terms of 
count. True peaks also have to be at least one pixel distant from 
a neighboring peak. Once the true “child” peaks are found over a 
detected region, a circular template is built for the child by using 
the smaller value of two pixels on the opposite sides of the peak. 
A weight is given to each child and a cost function is defined and 
minimized to deblend the sources (see Figure 1 of Lupton 2015, 
for a ID demonstration). 

In choosing a true peak, the area belonging to the peak and 
the count in it is ignored and a peak pixel is defined as true only 
if that pixel alone is sufficiently above its immediate valley or 
saddle point (compare to the discussion in Section 3.2.1). The 
very important technical issue of how to find the saddle point 
count is not discussed in that report. One drawback is that this 
selection criteria will miss real peaks that are below this large 
threshold but are not due to noise, for example with a very fiat 
profile. In this approach, deblending is completely independent 
of segmentation. Such pure deblending will also suffer when the 
detected regions cannot be fitted (the cost function minimized) 
with the combinations of templates over the detection. 

SExtractor uses a concept very similar to contour maps known 
as multi-thresholding. The count range of the detected region is 
divided into several layers. Then the separate objects in those 
layers are analyzed and “real” peaks are modeled to find the sep¬ 


arate local maxima and their contribution to the total count (see 
Figure 2 in Bertin and Arnouts 1996). clumpfind (Williams et al. 
1994) also uses a very similar approach. Therefore in this ap¬ 
proach, first a segmentation map is derived and used to model 
the objects, and then deblending is done. 

In SExtractor’s implementation, a fixed number of layers is 
defined for all objects, specified by deblend_nthresh in Section 
D. Therefore brighter objects, with a very large difference be¬ 
tween their peak count and the Sky, will receive layers that are 
much more widely separated than fainter objects. This will over¬ 
segment the fainter objects where the spacing between the layers 
is less and undersegment the brighter ones. The scale on which 
the count is divided will also significantly alter the result. In 
SExtractor, detection_type specifies the scale. If set to ccd, an 
exponential scale is used and setting it to photo will use a linear 
scale. Each will output significantly different results. 

Once peaks are found between the layers, they are accepted 
as a true detection or rejected as false based on the parameter 
DEBLEND_MiNC0NT which Specifies the fraction of the total inten¬ 
sity in that peak, or clump, compared to the total, undeblended, 
detection. Similar to the layer spacing problem mentioned above, 
this parameter is by definition, biased toward detecting more 
clumps in fainter objects. Since the threshold to accept a clump 
differs from one detection to another (based on the total count 
within each parent detection), a fixed peak will be detected in a 
fainter object and discarded in a bright one. Therefore when this 
selection criteria is adopted for a comparison of clumps in galax¬ 
ies, added with the layering bias mentioned above, it will bias the 
results by preferentially detecting more segments (or clumps) in 
objects with a lower total count before deblending. Another con- 
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Figure 18: SExtractor sensitivity test. Column 1: the input images, which are 
the same as those in Figure 12. Column 2: SExtractor’s segmentation maps. 
Column 3: SExtractor’s detections extended to three times the Kron radius. The 
SExtractor conhguration parameters can be seen in Appendix D. 

sequence of this method in finding true clumps is that the result 
will vary significantly depending on detect_thresh. When the 
threshold is high, the detected total count will be less, and there¬ 
fore more peaks will be detected in a given object. 

The highlighted (magnifying glasses and red arrows) regions 
in Figure 17 (for mock galaxies) and the large fraction of fiagged 
detections in Figure 19 (for real images) show how SExtractor’s 
deblending algorithm has failed in the low surface brightness re¬ 
gions. The primary reason for this failure is that SExtractor was 
designed for high thresholds. When reliable peaks are found 
(with a high threshold), each pixel is assigned to a peak, assum¬ 
ing all the peaks are a 2D Gaussian far larger than the region 
above the threshold. This is useful in cases where the peaks’ 
S/N and the threshold are sufficiently high and all peaks follow a 
Gaussian profile over the whole image. 

In practice, hardly any non-stellar objects of astronomical 
interest has a Gaussian profile and postage stamps are usually 
chosen to be larger than the objects in order to get an accurate 
estimate of the Sky value. Figure 17 shows two ways that SEx¬ 
tractor fails in deblending as the threshold decreases: (1) some 
of the very far detections of Figure 17(c.l) and (d.l), designated 
with arrows have been identified as belonging to the main object 
and (2) a connected region can be shared between two objects 
that are not physically connected to it as shown in the magnified 
regions of the same figure. These problems might not cause a 
significant problem when the magnitude of the very bright mock 



Figure 19: Tests of SExtractor on the same input images as Figure 13. First 
column: input image. Second column: SExtractor’s segmentation maps. Third 
column: 3,3.r^ elliptical apertures (see Section B.4). In the real images, because 
of their complexity, only the borders of the elliptical apertures are shown. In 
the mock image, the elliptical apertures are hlled. Most of the central object 
detections in a.3, b.3 and d.3 were flagged by SExtractor for possibly biased 
magnitude measurements. 

galaxy in this figure is desired, but when applied to low surface 
brightness and diffuse galaxies like Figure 19 (a)-(d), very seri¬ 
ous systematic errors can be produced. 

Figure 19 (a.l)-(d.l) shows SExtractor’s result for the same 
galaxies in Figure 13. Appendix D shows the configuration pa¬ 
rameters used to run SExtractor on these galaxies. The input PSF 
for these four real galaxies is obtained from the web interface of 
TinyTim^^ (Krist et al. 2011). The PSF of TinyTim is the theoret¬ 
ical PSF produced on the CCD. In order to create the processed 
images shown here, various steps of image processing have been 
applied to these images (see Koekemoer et al. 2007; Massey et 
al. 2010). Therefore, in practice, the actual PSF on the image 

http://tinytim.stsci.edu/cgi-bin/tinytimweb.cgi. On Chip 1, pixel position 
(2047, 1048), Alter F814W, Spectrum value 12, PSF diameter 3", focus 0.0. 
SExtractor does not accept kernels larger than 31x31 pixels. Therefore, the 
central region of this size was cropped from TinyTim’s output for input into 
SExtractor. 
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will be slightly different from that produced by TinyTim (see van 
der Wei et al. 2012). While this difference plays a very important 
role in model htting, for the demonstrative purposes on detection 
here, the TinyTim output is adequate. 

It is clear from Figure 19 that the two more diffuse objects 
in Figure 19(b) and (d) have been deblended into many more 
separate objects than the two brighter cases of Figure 19 (a) and 
(c). As explained above, this is due to the inherent biases of 
SExtractor’s deblending algorithm. 

B.3. Combining Detection and Segmentation 

A discrete hierarchical Markov image modeling approach 
was introduced in astronomy recently by Vohmer et al. (2013). 
This is a fundamentally different approach to detection and seg¬ 
mentation. It is reviewed here because like the more popular 
methods used today, the detected region needs to be grown to 
complete the detection as we will discuss next in Section B.4. In 
this approach object detection and segmentation are done simul¬ 
taneously. The basic technique comes from image processing 
applications in particular the quadtree of Laferte et al. (2000). 

Hierarchical Markov models have the advantage that they al¬ 
low to use multi-band images at potentially different pixel reso¬ 
lutions without the need to modify the data as most techniques 
used in astronomy currently need to do. This makes them “scale 
free.” However Vohmer et al. (2013) acknowledge that their im¬ 
plementation (MARSIAA) does not account for PSF variation be¬ 
tween the inputs. In this approach the problem is labeling of 
the image pixels “which are spatially connected through a pre- 
dehned neighborhood system. These labels correspond to dis¬ 
crete classes of objects with similar surface brightness”(Vohmer 
et al. 2013, page 3). Using the multi-band data as input, the pro¬ 
cedure iterates along the hierarchies to hnd a distribution for the 
labels over the image, see Laferte et al. (2000) for details. 

In this approach there is no specihc functional basis to do the 
detection, therefore it is considered as a model-independent de¬ 
tection and segmentation technique. However, instead of known 
functions, it classihes the image pixels into “classes” which are 
shared by ah the objects in the image. For example, noise has a 
class of 0 and as the gradient over the object increases, the class 
also increase; see Figure 9 in Vohmer et al. (2013) for how the 
class changes toward the center of the object. The number of 
classes is a free parameter specihed by the user. The hnal simul¬ 
taneous detection/segmentation depends on the total number of 
classes and the dynamic range of the image. To detect low sur¬ 
face brightness galaxies, Vohmer et al. (2013) had to truncate or 
clip the image pixels beyond 20a, therefore this approach cannot 
simultaneously detect/segment very bright and very faint objects 
in one image. 

An interesting consequence of this technique is that it does 
not need to dehne a count threshold, however, it can only get 
^ 1.5a close to the Sky value. Figures 8-10 of Vohmer et al. 
(2013) show the results of this technique compared to SExtrac¬ 
tor’s segmentation maps for three galaxies. In the three cases, the 
detected area of SExtractor’s segmentation map with a threshold 
at 1.5 a was comparable or even larger (better covering the ob¬ 
jects). This technique heavily relies on iterative processes and 
thus can be very slow. Vohmer et al. (2013) report it was at least 
90 times slower than SExtractor for the same data set. 


B.4. Growing Detections Below the Threshold 

The high thresholds that were used to avoid false detections 
cause a very large fraction of the object to not be initially de¬ 
tected (see Figure 17(a.l)-(d.l)). Failing to detect a bright ob¬ 
ject’s complete area will deprive us of ah the valuable infor¬ 
mation lying hidden there and cause a systematic under- (over- 
)estimation of the object’s magnitude (Sky value). Therefore 
when using such large thresholds, it is necessary to use the re¬ 
gions above the threshold as a basis to “grow” the detections and 
include their fainter parts. 

The Kron radius, (Kron 1980), is one of the most popular 
methods to grow the detection into regions of the image with a 
count lower than the threshold. It is dehned as the light-weighted 
average radius of a prohle. 


From Figure 17 it is clear that the elliptical parameters and 
rk are found based on the region above the threshold, therefore 
the fainter regions of the galaxies that often harbor valuable in¬ 
formation on the dynamical history of the galaxy and might have 
different morphological parameters cannot be directly measured 
in the existing signal based approach to detection. 

It is generally considered that contains > 90% of the to¬ 
tal count, integrated to inhnity (F}), of a galaxy light prohle (see 
Kron 1980; Infante 1987). Therefore in the literature it is com¬ 
mon to use 2.5rk to hnd the total count of a galaxy. In order 
to test this claim, SExtractor’s flux_auto is divided by for 
the four different masked 3,3.ry^ apertures shown in Figure 17. 
For this prohle, 3,3.ry^ contains 74.30%, 78.21%, 80.22% and 
87.33% of the total count. The area within the 3,3.has been 
masked to show the extent of the object that has not been de¬ 
tected. The elliptical apertures containing 90% and 95% of F} 
are also shown in Figure 17 for visual comparison. 

Graham and Driver (2005) did a thorough analytic study of 
the mathematical (continuous) ID Sersic prohle to show how 
such differences from the expected > 90% could occur in the 
calculation of the Kron radius. They showed how dependent the 
total count within the 2.5 times the Kron radius is to the area that 
is used to hnd the Kron radius. If this area is extended to inhnity, 
it is indeed as expected and 2.5 times the Kron radius contains 
> 90% of F/, even for ^ = 10 prohles. However, integration 
to inhnity is not a realistic condition. In practice any measure¬ 
ment of the Kron radius is conhned to a certain aperture. By 
considering the effects of the hnite aperture that is used in real 
measurements of Vk, they show how the values reported above, 
which have also been reported by other authors are reasonable 
(see Bernstein et al. 2002; Benitez et al. 2004). 

This discussion demonstrates that the statement that an aper¬ 
ture > 2rk will contain nearly ah (> 90%) of the light and area 
of a galaxy is not generic and is only practically correct for very 
hat prohles. If ignored and used for ah galaxies in a survey or 
image, it can result in a very signihcant underestimation of the 
total light and boundaries of the sharper objects, even if it is 
a pure ellipse as in Figure 17. The underestimation will cause 
an overestimation in the Sky value of the image, further com¬ 
pounding the resulting systematic bias. The mock galaxy used 

Fj is approximated by the total count within 30rg. Note that the larger ellipse 

in Figure 17 which contains 95% of the total light is on Vr^. 
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in this example has a Sersic index {n) of 4.00, while in the local 
universe, galaxies with n as large as 11.84 have been reported 
(Kormendy et al. 2009). The case for distant galaxies becomes 
much worse, since their apparent in pixels, due to the angu¬ 
lar diameter distance, is much smaller with much lower surface 
brightness. 

Figure 18 shows SExtractor’s output for mock galaxies, both 
flat, low Sersic index, and sharp or high Sersic index. Even with 
DETECT_THRESH=l, only a few pixels have been found for each 
object, the elliptical parameters, namely center, axis ratio, posi¬ 
tion angle and major axis length {vk) are all measured using these 
few pixels. The second argument to phot_autoparams is set to a 
minimum radius such that if Vk is too small, this minimum radius 
is used instead of it. This is why most of the ellipses in Figure 18 
have similar sizes. However, such hand input values could add 
more bias to the flnal result (see Section 5.3 for a discussion on 
the dispersion of the magnitudes measured in this method). 

By deflnition, depends on flnding the full ID light profile 
of an object. In order to do that, it is necessary that the cen¬ 
ter of the profile, its position angle and axis ratio are accurately 
defined. Furthermore, a good enough sampling of the profile at 
various (r, 0) is needed in order to approximate the count at a 
certain radius. These preconditions are only valid in the most 
optimistic cases, like the bright purely elliptical mock profile 
of Figure 1(b). In general only a very rare fraction of galaxies 
display such a simple and bright profile, clearly separated from 
neighbors or image edges. 

Galaxies can have multiple components, for example, the 
bulge, a bar, and a disk. Each can have a different axis ratio and 
position angle. The central parts (bulge and a possible bar) are 
much brighter than the outer disk. Therefore their effect on the 
overall, light-weighted, estimation of ellipticity will completely 
outweigh the fainter outer parts, which might not even be de¬ 
tected because of the high thresholds (detect_thresh>i) that are 
commonly used. This compounds the problem mentioned in the 
previous paragraph for a real galaxy’s total light because all the 
elliptical parameters in one profile can change in a real galaxy 
due to their rich dynamical history. 

The four real galaxies in Figure 19 visually appear to have 
a distinct region in the noise. Whether they are composed of 
separate objects in the same line of sight or are actually one 
object cannot be judged from one image. In all four exam¬ 
ples, because the connectivity of the objects is below the thresh¬ 
old (detect_thresh=i), SExtractor considers them to be spatially 
separated objects. Furthermore, in none of these examples can 
the connectivity be modeled as a simple ellipse to be able to 
model. 

Figure 19(a.l) appears to have one bright clump and a dif¬ 
fuse host. Because of the relatively flat, diffuse structure, the 
ellipses showing 3,3.ry^ for the three detections have become ex¬ 
tremely large. When this happens, the flnal result is dependent 
on the accuracy of Sky subtraction, because a large area of sky is 
also included as part of the object. The large fraction of flagged 
detections in the real galaxies of Figure 19 shows how mag_auto 
has been flagged as unreliable for most of their sub-components. 
A user of SExtractor will thus have to either ignore flagged detec¬ 
tions or rely on mag_iso, using only the segmented regions above 
the threshold. If they choose the latter, for such diffuse objects, 
MAG_iso is going to miss a very large fraction of the object’s area 
and photon count. 

Due to the sheer number of galaxies that researchers use for 


their studies, it is not feasible to check each detection by eye 
to check when mag_auto should be used or mag_iso. SExtractor 
does provide mag_best to automatically choose between the two, 
but it is not commonly used in studies today. The main reason is 
that MAG_AUT0 and mag_iso are not consistent and cannot simply 
be compared with each other in the study of a large number of 
galaxies. 

In Figure 19(e), because of the strong gradient created by 
bicubic-spline interpolation in SExtractor’s Sky measurement (see 
Figure 16(1)) no pixel from the the profile on the corner of the im¬ 
age has been detected. Looking at SExtractor’s detections in Fig¬ 
ure 19 (e.3), it is clear that for such objects, lying on the edges of 
images, any modeling method of trying to account for their faint 
regions fails. A very large portion of their total count is left out of 
their designated Kron aperture. Like the case for clumpy objects, 
all the necessary requirements to define fail when enough of 
the object is beyond the edge of an image. Such objects have to 
be taken into account, particularly if they are the wings of bright 
objects. While the total magnitude of such objects would always 
be flagged by software like SExtractor, and not included in any 
scientific work, failing to account for as much of their light as 
possible will systematically increase the Sky value that is calcu¬ 
lated in such an image. Such cases are common when postage 
stamp images (for example. Figure 19 (a.l)-(d.l)) of a selected 
group of galaxies in a larger survey are studied. 

The discussions in this section show that the only way to get 
as close as possible to the “true” total photon count and area 
of an object of interest is to use lower thresholds. Trying to 
correct for the photon count which is lost, when a high threshold 
that is required in the signal-based detection technique is used, 
will only be a source of systematic bias in measurements and 
the resulting scientific results. The problems mentioned in this 
section were the primary motivation in flnding a new approach 
to detection and segmentation that would mitigate many of the 
problematic issues described above. 

C. FINDING THE MODE 

In Section A.l, the existing methods to And the mode of an 
image were reviewed. Here a novel approach to flnding the mode 
of a distribution is proposed. It is very accurate as long as the 
mode of the distribution is approximately symmetric around it, 
like the pixel count distributions of Figure 1 . This is a safe con¬ 
straint for images of astronomical targets, which generally have 
more fainter pixels than brighter ones (the yellow histograms of 
Figure 1). When the noise is significant, the noise will ensure 
a symmetric mode similar to the examples of Figure 1 . As the 
signal becomes more significant, the quantile of the mode will 
shift to lower values. 

Given an ordered data set X with n elements, where Xn > 
Xn-i{n > 1). A mirrored distribution is defined about the mirror 
point located at index m as follows. All the data prior to m are 
placed in opposite order after m, as if a mirror was placed there. 
The elements of the mirrored data set are thus defined with 

Xi i < m 

Xjfi + {Xjfi i ^ m 

The mirrored distribution is therefore identical to the orig¬ 
inal distribution until and including Xf^. M has 2(m — 1) + 1 
elements and like X, it is ordered. Figure 20 shows one example 
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Figure 20: Mirrored distribution compared with the original distribution. The data set is the pixels of Figure 9(b.2) multiplied by 10^ seconds. Mm is the value on 
the mirror point. Blue and green histograms are the histograms of the original and mirrored distributions, respectively. The thick and thin lines are the cumulative 
frequency plots of the original and mirrored distributions, respectively. The percentage value shown after the label of each plot shows the percentile of the original 
distribution on which the mirror was placed. It is clear that as the mirror point gets closer to the mode of the distribution, the two lines diverge more smoothly. The 
cumulative frequency curve is scaled to the tip of the blue histogram. 


data set (the pixels of Figure 9(b.2)) with mirrored distributions 
overplotted. The mirror points in Figure 20 are placed at various 
quantiles of the original distribution in each case. The histogram 
is only displayed in this demonstration to aid the understanding. 
It is not used at all during the analysis. 

Figure 20 shows that when the mirror is very far from the 
mode, the cumulative frequency plot of the mirrored distribu¬ 
tion will deviate below or above that of the original pixels very 
sharply. As the mirror approaches the mode, the maximum dif¬ 
ference of the two cumulative frequency plots in the vicinity of 
the mode becomes smaller. Therefore, on the symmetric mode 
of a distribution, the two cumulative frequency plots should re¬ 
main very similar for a large number of pixels (or distance in the 
ordered array) after the mirror. 

Let 5(M/,m) represent the difference between the cumula¬ 
tive frequency plots of the original and mirrored data sets for M/, 
when the mirror is placed on the index m. For any given mir¬ 
ror position and positive integer A, let A(m) be the maximum 
|5(M/,m)| for all m < / < m + A. Then the symmetric mode of 
a distribution can be found by finding m such that A(m) is min¬ 
imized. The Golden section search algorithm (Kiefer 1953) is 
used to find the mirror point that minimizes A(m). 

If the distribution has several symmetric local maximums, 
and as long as A is large enough, the first (with the lowest quan¬ 
tile) will be found, even if the brighter peak has more pixels in 
its proximity (see Figure 21(c) for a distribution with two local 
maximums in the histogram). This is exactly what is desired in 
astronomical image processing because if there is a second local 
maximum with a larger count and larger number of pixels, then 
that local maximum is not due to noise. Noise will always pro¬ 
duce the first (lowest count) local maximum, which can be called 
the mode. In most astronomical applications of low S/N objects, 
the first local maximum will be the only one. If not, it will have 
more pixels in its vicinity than any other local maximum in the 
pixel distribution. 

The mirror cumulative frequency plot found for the mode 
should be approximately equal to the actual data in the vicinity 
of the mode. Beyond that, it should ideally always remain below 
the actual data (see Figure 21). However, in minimizing A(m) 
as defined above, the resulting mirror cumulative frequency plot 
will zig-zag around the actual distribution since |5(Mi,m)| was 
used. The cumulative frequency plot is fundamentally a count¬ 
ing function. Therefore, the error in choosing an index, say m, 
follows a Poisson distribution. Thus the standard deviation in 
finding m is Hence a limit can be set on how negative 

can become. Taking a to be a given constant speci¬ 
fied by the user, if < —a^/m, A(m) is set to a fixed 


constant, let’s call it The golden section search is modified 
from its standard algorithm to choose the interval with the low¬ 
est values when it confronts C. Recall that such conditions only 
occur when m is a significant overestimation of the mode. 

Figure 21(a) shows the result for the same distribution of 
Figure 20. In Figure 21(b) this technique is applied to the more 
skewed distribution of the convolved image of Figure lO(c.l). 
A very bright spiral galaxy has caused the postage stamp pixel 
distribution to be more skewed. Figures 21(c) and (d) show the 
distributions of two mock images with two very bright and large 
Sersic profiles of ^ = 2.5 and n = 0.5 placed in a small postage 
stamp respectively. These are two extreme cases, showing how 
the quantile of the mode can reach extremely low quantiles. 

As the distribution becomes more skewed, the mode shifts 
to lower quantiles; see Figure 21. This is the basis of the new 
method to define the threshold of the image (explained fully in 
Section 3.1.2 and Section 4.2.1). In these examples, the bound¬ 
aries for the golden section search are set to the 0.01 and 0.51 
quantiles of the image. As long as data is present in the image 
(which create a positive skew) the mode will not exceed the me¬ 
dian, thus the higher limit. In these examples, A is set to be 1/2 
of the total number of pixels and in order to be more efficient, 
5(M/,m) is sampled in 1000 equally spaced points between Mm 
and Mm+N- oc is also taken to be 1.5. 

Depending on the brightness of the data pixels and the frac¬ 
tion of pixels that they occupy, the mode calculated with this 
technique can loose its accuracy or it can cease to exist. Also, 
in extremely skewed distributions, the golden section intervals 
might miss the mode and converge on a lower mirror point. There¬ 
fore a measure of quality is necessary to accompany the mode. 

Let a be the value of the 0.01 quantile of the mirrored dis¬ 
tribution. The minimum is not chosen due to its scatter. Take 
h to be the value on the mode and finally, let c be the value of 
the first Mi beyond the mode where |5(M; ,m) | > a ^/m, same a 
as above. This point is marked on the plots of Figure 21 with a 
dotted vertical line. These three values can define a measure of 
“symmetricity” (s) about the mode with: 

c — b 


The symmetricity for the examples of Figure 21 are respec¬ 
tively: 0.43,0.25,0.84 and 0.74. Based on these and other simu¬ 
lations, 5- > 0.2 appears to be a good symmetricity of an accurate 

The value of C is an internal issue and is irrelevant to this discussion or the 
user. It is only necessary that it should be larger than the largest possible 
A(m). 
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Figure 21: Finding the mode by minimizing A(m). Similar to Figure 20, with the mirror placed on the mode. The dotted line shows the first value where |6(M/,m) | > 
1.5 ^/m, which is used to find the “symmetricity” of the mode (see text), (a) Same data set as Figure 20. (b) The convolved image of Figure lO(c.l): a bright and large 
galaxy skewing the distribution, further shifting the mode to lower percentiles, (c) The convolved and noisy image of a circular —20 magnitude mock Sersic profile 
with n — 2.5, = 50 pixels, truncated at Ivg with center on (100, 140) in a 200 x 200 pixel image, (d) Similar to (c) but —28 magnitude with n = 0.5, = 100, 

truncated at Sve placed at (110, 285). 


“symmetric” mode, when the mode is above the ^ 0.02 quantile 
of the image. When the output mode approaches such quantiles, 
b —am the equation above will be too small to give an accurate 
result. 

In addition to the initial cost of sorting of the data set, the 
rest of the process is very simple and not CPU intensive, since 
only 1000 points are checked in each round of the golden sec¬ 
tion search. Therefore the CPU cost of running this technique is 
similar with the cost for a-clipping which also requires sorting, 
to find the median, but needs multiple passes through the whole 
data. 

D. SEXTRACTOR PARAMETERS 

Below is a list of the S Extractor parameters used in this study. 
We make no claim that these are “the best” (or worst) possible 
set of parameters. The values for each parameter were selected 
based on values reported by other research teams which exten¬ 
sively used SExtractor in their astrophysical data analysis and its 
manual; see Appendices A and B. 

The file PSF. conv refers to the PSF used to create mock im¬ 
ages when SExtractor was applied to the mock images and to 
that produced by TinyTim when real HST images were used; see 
Appendix B for more details. 


MEM0RY_0BJSTACK 

3000 

MEM0RY_PIXSTACK 

300000 

MEM0RY_BUFSIZE 

1024 

DETECT_TYPE 

CCD 

GAIN_KEY 

GAIN 

PIXEL_SCALE 

1.0 

FITS_UNSIGNED 

N 

BACK_TYPE 

AUTO 

BACK_SIZE 

64 

BACK_FILTERSIZE 

3 

BACKPH0T0_TYPE 

LOCAL 

BACKPH0T0_THICK 

24 

WEIGHT_TYPE 

BACKGROUND 

WEIGHT_GAIN 

N 

FILTER 

Y 

FILTER_NAME 

PSF.conv 

MASK_TYPE 

CORRECT 

INTERP_MAXXLAG 

16 

INTERP_MAXYLAG 

16 

INTERP_TYPE 

ALL 


DETECT_THRESH 

1 

ANALYSIS_THRESH 

1 

THRESH_TYPE 

RELATIVE 

DETECT_MINAREA 

5 

DEBLEND_NTHRESH 

64 

DEBLEND_MINCONT 

0.001 

CLEAN 

Y 

CLEAN_PARAM 

1.5 

FLAG_IMAGE 

flag.fits 

FLAG_TYPE 

OR 

SEEING_FWHM 

1.2 

starnnw_name 

default.nnw 

MAG_ZER0P0INT 

0.0 

PHOT_APERTURES 

5 

PH0T_AUT0PARAMS 

3,3.5 

PH0T_AUT0APERS 

0,0 

PHOT_FLUXFRAC 

0.5 

PHOT_PETROPARAMS 

2.0,3.5 

SATUR_LEVEL 

50000.0 

SATUR_KEY 

SATURATE 

CHECKIMAGE_TYPE 

SEGMENTATION,APERTURE S,BACKGROUND 

CHECKIMAGE_NAME 

seg.fits,aper.fits,back.fits 

CATALOG_NAME 

SEresults.txt 

CATALOG_TYPE 

ASCII_HEAD 

PARAMETERS_NAME 

./parameters/sextractor/sextractor.param 

VERBOSE_TYPE 

NORMAL 

E. NOISECHISEL PARAMETERS 

The full list of parameters that was used for Noise Chisel in 
this paper is printed here. If a different parameter was used, it is 
explained in the text. Note that we do not consider these param¬ 
eters to be “the best” for any generic data set. Eor the processing 

of Eigure 11, — 

nchl=4, —lmesh=256 and the following two 

on/off options were also activated — fullinterpolation and 

—fullsmooth. 

activated. 

Eor the HST images, — skysubtracted was 

# Input: 


hdu 

0 

khdu 

0 

skysubtracted 

0 

minbfrac 

0.5 
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minnumfalse 

100 

# Mesh grid: 

smeshsize 

32 

Imeshsize 

200 

nchl 

1 

nch2 

1 

lastmeshfrac 

0.51 

mirrordist 

1.5 

minmodeq 

0.49 

numnearest 

10 

smoothwidth 

3 

fullconvolution 

0 

fullinterpolation 

0 

fullsmooth 

0 

# Detection: 

qthresh 

0.3 

erode 

2 

erodengb 

4 

opening 

1 

openingngb 

8 

sigclipmultip 

3 

sigcliptolerance 

0.2 

dthresh 

-0.1 

detsnminarea 

15 

detsnhistnbins 

0 

detquant 

0.99 

dilate 

3 

# Segmentation: 

segsnminarea 

25 

segquant 

0.99 

segsnhistnbins 

0 

gthresh 

0.5 

minriverlength 

15 

objbordersn 

1 
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